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Summary 


Details  are  given  of  an  implicit  six-point  finite-difference  scheme 
for  solving  two- temperature  chemical  nonequilibrium  laminar  boundary- layer 
flows  in  ionizing  argon.  The  analysis  extends  previous  work  by  considering  the 
radiation- energy  loss  and  the  chemical  reactions  in  the  plasma  of  the  ionizing 
boundary  layer.  The  variations  of  transport  properties  based  on  the  known 
elastic- scattering  cross-sections  for  an  argon  plasma  across  the  boundary  layer 
are  considered.  The  effects  of  the  chemical  reactions,  radiation-energy  loss 
and  the  electric  sheath  on  the  boundary- layer  structures  are  discussed.  Both 
the  flat-plate  and  the  shock- tube  sidewall  boundary  layer  flows  are  analyzed 
and  compared  with  interferometric  data  obtained  using  the  UTIAS  10  cm  x 18  cm 
Hypervelocity  Shock  Tube  at  shock  Mach  numbers  Mg  ^ 13  and  ~ l6  at  an  initial 
argon  pressure  pQ  ~ 5 torr  and  temperature  T0  ~ 300  K.  Fairly  good  agreement 
was  obtained  between  analysis  and  experiment  for  both  types  of  boundary  layers. 
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1.  INTRODUCTION 


1.1  General  Considerations 

An  understanding  of  boundary  layer  flows  in  partially  ionized  gases 
is  helpful  in  designing  successful  spacecraft  for  re-entry  into  the  Earth' s 
atmosphere  at  hypersonic  conditions.  It  also  provides  insight  into  the  physical 
phenomena  of  interactions  between  solid  surfaces  and  plasma  flows.  The  presence 
of  ions  and  electrons  introduces  new  .transport  mechanisms  and  chemical  reactions 
in  the  boundary  layer.  The  magnitude  of  the  various  transport  properties  of  an 
ionized  gas  can  be  markedly  different  frcm  a perfect  gas.  Therefore,  the  boundary 
layers  in  an  ionizing  gas  are  generally  more  complex  than  those  in  noniomzed-gas 
flows.  Despite  years  of  research,  .boundary-layer  flows  of  partially-ionized  gases 
are  not  fully  understood  experimentally  and  theoretically. 

The  character  of  the  ioniziig  boundary  layer  problem  was  schematically 
described  by  KnOOs  (Ref.  l) . The  following  characteristics  are  important  in 
considering  partially-ionized  boundary-layer  flows: 

(a)  Transport  properties 

(b)  Interactions  between  moving  plasma  and  metal  surface 

( c ) Atomic -collision  processes 

(d)  Chemical  reactions 

(e)  Radiation-energy  transfer 

(f)  Electromagnetic  fields 

The  full  boundary  layer  problem  is  exceedingly  complex,  and  only  a 
few  cases  have  been  treated  by  early  investigators.  Usually,  same  approximations 
are  made  to  suit  a given  problem  and  to  reduce  the  computation- time  costs. 

In  general,  a mixture  of  an  ionizing  gas  is  composed  of  molecules, 
atoms,  molecular  ions,  atomic  ions  and  electrons.  However,  since  the  dissocia- 
tion energy  is  much  less  than  the  ionization  energy,  ionization  can  be  considered 
to  become  appreciable  only  after  dissociation  is  practically  completed.  There- 
fore, the  mixture  is  assumed  to  be  composed  only  of  atoms,  atomic  ions  and 
electrons.  The  presence  of  electrons  in  a gas  introduces  some  features  quite 
different  from  those  encountered  in  chemical  dissociations.  For  example,  the 
collisional  energy- transfer  processes  between  electrons  and  heavy  particles 
(atoms  and  ions)  are  relatively  slow,  giving  rise  to  the  possible  situation 
that  the  electrons  may  have  a temperature  much  different  from  that  of  the  heavy 
particles.  The  extremely  low  mass  of  the  electrons  yields  a species  possessing 
a thermal  conductivity  that  can.be  much  greater  than  that  of  the  other  species. 

When  such  a gas  is  in  contact  with  a cold  surface , a space- charge  sheath  is  formed 
which  may  affect  the  energy  transfer  to  the  surface.  The  electrons  may  have  a 
higher  temperature  than  the  heavy  particles  near  the  cool  surfaces.  In  such 
cases,  the  electrons  make  a greater  contribution  to  the  electrical  and  thermal 
conductivity  than  would  be  expected  solely  on  the  basis  of  their  numiber  density. 

Finally,  the  charged  species  are  sensitive  to  electromagnetic  fields  yielding 
a possible  method  of  controlling  the  energy  transfer  processes  between  electrons 
and  ions.  Therefore,  the  boundary- layer  flows  in  ionizing  gases  are  exceedingly 
more  complex  than  in  nonionized  or  dissociated  gases. 

The  state  of  the  mixture  of  atoms,  ions  and  electrons  is  uniquely 
described  by  three  independent  state  parameters:  pressure,  temperature  and 
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degree  of  ionization.  In  general,  the  mixture  of  an  ionizing  gas  is  in  a state 
of  nonequilibrium,  that  is,  both  thermal  nonequilibrium  and  chemical  nonequili- 
brium. In  thermal  equilibrium,  which  in  general  cannot  be  expected  to  occur  in 
a boundary-layer  flow,  the  temperatures  (Ta  and.  Te)  of  both  heavy  particles  and 
electrons  are  equal.  In  chemical  equilibrium,  the  degree  of  ionization  a is 
immediately  adjusted  to  its  local  equilibrium  value,  so  that  the  degree  of 
ionization  can  be  given,  as  a function  of  pressure  p and  temperature  T via  the 
Saha  equation.  In  the  frozen  state,  the  electron  number  density  production 
rate  ne  is  equal  to  zero.  The  following  models  have  been  considered  by  a 
number  of  authors  in  solving  ionizing  boundary  layer  flows: 


One-temperature  equilibrium: 
Two- temperature  equilibrium: 
One- temperature  frozen: 

Two- temperature  frozen: 


Te  = Ta}  a = f(p>  Ta) 

Te  t Ta»  a = f(p’  Te  or  Ta^ 

Te  = Ta,  ne  - 0 

T / T , n =0 


e ' a"  e 

One-temperature  nonequilibrium:  T = T , n ^ 0 

6 a 6 

Two-temperature  nonequilibrium:  T ^ T , n ^ 0 

g a g 

The  thermodynamic  quantities  and  the  descriptions  of  equilibrium, 
frozen  and  nonequilibrium  flows  are  given  in  Appendix  B. 


The  following  brief  review  may  be  helpful . Many  investigators  have 
treated  weakly  ionized,  colli si on- dominated  boundary  layers.  Their  main  aim 
was  to  study  the  effects  produced  on  the  electrical  characteristics  of  Langmuir 
probes.  Examples  are  the  incompressible  flow  of  a weakly  ionized  gas  treated 
by  Su  and  Lamb  (Ref.  2)  and  the  Couette-flow  problem  studied  by  Chung  (Ref.  3). 
The  kinetic  theory  of  ionized-gas  flows  was  used  in  the  analysis.  Recently, 
Chung,  Talbot  and  Touryan  (Ref.  4)  have  summarized  the  theoretical  results  for 
electric  probes. 


Based  on  thermal  equilibrium  in  temperature  and  chemical  reactions, 

Fay  and  Kemp  (Ref.  5)  have  studied  the  heat  transfer  to  a shock- tube  end-wall 
from  an  ionized  monatomic  gas  and  Knoos  (Ref.  l)  generalized  it  to  a simple 
thermal  Rayleigh  boundary  layer  in  an  equilibrium  flow.  Back  (Ref.  6)  studied 
the  heat  transfer  through  a one-temperature  laminar  boundary  layer  from  a par- 
tially-ionized gas  to  a highly- cooled  wall  for  frozen  and  equilibrium-flow 
models  based  on  similar  assumptions . A finite- difference  method  was  applied 
by  Blottner  (Ref.  7)  to  a one-temperature  nonequilibrium  laminar  boundary  layer 
in  ionizfag  air.  Park  (Ref.  8)  analyzed  the  frozen  and  equilibrium  flow  over  a 
flat  plate  and  at  an  axi symmetric  stagnation  point  based  on  similar  and  one- 
tenperature  models.  Fins on  and  Kemp  (Ref.  9)  extended  the  theory  of  Fay  and 
Kemp  to  stagnation-point  heat  transfer.  Using  one- temperature  and  constant 
transport  properties,  the  equilibrium,  frozen  and  nonequilibrium  solutions  were 
obtained  by  Liu  (Ref.  10)  through  an  integral  method. 

For  the  two-temperature  boundary  layer,  Sherman  and  Reshotko  (Ref.  11) 
have  obtained  the  electron  temperature  profiles  for  chemical- equilibrium  flow 
based  on  similar  solutions.  Nishida  and  Matsuoka  (Ref.  12)  solved  the  similarity 
equations  for  frozen  flow  with  constant  transport  properties.  Analyses  of  flat- 
plate  boundary  layers  in  parti ally- ionized  gases  with  thermal  nonequilibrium  and 
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recombination  were  investigated  by  Tseng  and  Talbot  (Ref.  13)  based  on  similar 
solutions  and  constant  transport  properties.  Recently,  Takano  and  Akamatsu 
(Ref.  l4)  used  a finite- difference  method  to  solve  the  shock- tube  side-wall 
boundary-layer  flow  with  constant  transport  properties.  The  thermal  Rayleigh 
boundary  layer  flow  for  partially-ionized  argon  with  varied  transport  proper- 
ties were  studied  numerically  and  experimentally  by  Mansfeld  (Ref.  15)  for 
thermal  and  chemical  nonequilibrium  cases.  Honma  and  Karauro  (Ref.  16)  studied 
an  ionizing  nonequilibrium  boundary  layer  behind  a moving  shock  wave  by  using 
a finite-difference  scheme. 

The  numerical  methods  for  solving  the  boundary  layer  equations  can 
be  divided  into  the  following  categories: 

(a)  Local- similarity  method 

(b)  Integral  method 

(c)  Difference-differential  procedure 

(d)  Series-expansion  method 
. (e)  Perturbation  method 

(f)  Finite-difference  method 

With  the  exception  of  the  finite- difference  scheme,  all  these  tech- 
niques involve,  in. one  way  or  another,  the  reduction  of  the  nonlinear  partial- 
differential  equations  to  ordinary-differential  equations.  In  the  local- 
similarity  method  the  history  of  the  flow  is  ignored  except  insofar  as  it 
appears  in  the  calculation  of  the  variable  x (or  |),  where  x is  the  coordinate 
along  the  surface.  This  results  in  a set  of  ordinary-differential  equations 
with  two-point  boundary  conditions.  In  the  integral  method  one  or  more  assump- 
tions are  made  regarding  the  profiles  of  the  flow  quantities.  The  equations 
used  are  obtained  by  taking  suitable  integrals  of  the  boundary-layer  equations 
across  the  boundary  layer.  The  boundary- layer  equations  reduce  to  a system  of 
ordinary-differential  equations  of  the  initial-value  type.  In  the  difference- 
differential  procedure,  the  derivatives  in  the  direction  along  the  surface 
are  replaced  with  finite-difference  relations  and  the  nonlinear  partial- differ- 
ential equations  reduce  to  ordinary- differential  equations  with  two-point 
boundary  conditions.  In  the  series-expansion  method,  the  coefficients  of  a 
series  in  an  x-dependent  variable  are  obtained  from  a solution  of  ordinary- 
differential  equations.  The  expansion  variable  depends  on  the  external-flow 
conditions.  The  perturbation  method  is  based  on  the  concept  that  a perturba- 
tion of  a known  boundary- layer  solution  is  considered  and  an  expansion  is 
carried  out  in  terms  of  a parameter.  A critical  review  of  the  early  work  up 
to  1969  was  given  by  Blottner  (Ref.  17). 

Two  difficulties  exist  in  the  analysis  of  ionizing  boundary— layer 
flows:  (l)  the  evaluation  of  the  reaction-rate  coefficients  near  the  wall, 

(2)  the  boundary  conditions  for  the  degree  of  ionization  and  the  electron 
temperature  at  the  wall.  First,  near  the  wall,  where  the  temperature  of  the 
heavy  particles  is  in  equilibrium  with  the  wall  temperature , the  temperature 
of  the  heavy  particles  is  very  small  compared  with  that  at  the  edge  of  the 
boundary  layer.  Near  the  wall  the  electron-number  density  is  also  very  low. 

In  this  low  temperature  and  low  electron-number- density  domain,  thermal 
ionization  hardly  occurs.  Consequently,  thermal  transport  processes  will 
dominate.  The  reverse  chemical-reaction-rate  coefficients  for  atom-ion- 
electron  and  electron-ion-electron  collisions  are  extremely  large  and  difficult 
to  evaluate  in  that  domain.  Second,  the  boundary  conditions  for  the  degree  of 
ionization  and  the  electron  temperature  at  the  wall  are  usually  determined  from 
the  collision-free  sheath  theory.  However,  some  authors,  for  example  Mansfeld 
(Ref.  15),  found  that  owing  to  the  assumptions  and  incomplete  description  of  the 
electric  sheath  a comparison  of  theoretical  and  experimental  results  would  be 
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of  questionable  value.  Careful  experiments  for  electron-nunber  density  near 
the  wall  must  be  done  in  order  to  check  the  validity  of  the  electric- sheath 
theory. 


The  difficulty  in  using  the  finite-difference  method  for  an  ionizing 
boundary  layer  lies  in  the  stability  of  the  scheme  and  in  significant  computa- 
tion times.  The  stability  criterion  for  the  set  of  strongly-coupled  nonlinear 
partial— differential  equations  with  their  boundary  conditions  of  a mixed 
Neumann/Dirichlet  type  in  the  finite- difference  scheme  is  difficult  to  evaluate. 

In  order  to  avoid  the  difficulty  of  stability,  Mansfeld  (Ref.  15)  applied  the 
backward  implicit  method  in  the  time- dependent  one- dimensional  Rayleigh  problem. 
However,  his  program  is  near  the  maximum  acceptable  computation  time.  In  the 
present  two-dimensional  boundary-layer  flow,  which  is  more  complex  than  the 
Rayleigh  boundary-layer  flow,  the  stability  criterion  and  computation  time  should 
be  examined  carefully. 

Blottner  (Ref.  17)  mentioned  that  the  iteration  procedure  for  controlling 
the  nonlinear  terms  is  not  required  for  a dissociated  boundary- layer  flow.  However, 
when  the  variations  of  the  transport  properties  across  the  boundary  layer  are 
taken  into  account  in  the  ionizing  boundary- layer  equations,  a successive  iteration 
procedure  is  necessary  in  the  present  problem.  This  iteration  scheme  increases 
the  computation  time.  Therefore,  in  the  present  calculation,  the  implicit  six- 
point  finite- difference  method  and  nonequidistant  step  sizes  are  applied  in  order 
to  decrease  the  computation  time. 

1.2  Scope  of  Present  Work 

The  present  paper  gives  the  details  of  an  implicit  six-point  finite- 
difference  scheme  for  solving  the  nonlinear  partial— differential  equations  of 
thermal  and  chemical— nonequilibrium  boundary-layer  flows  in  ionizing  argon.  The 
transport  properties  evaluated  from  known  elastic- scattering  cross-sections  of 
the  plasma  are  varied  across  the  boundary  layer.  The  radiation- energy  loss  of 
the  plasma  and  the  appropriate  chemical  reactions  are  both  considered.  The 
flat-plate  and  shock- tube  sidewall  boundary- layer  flows  are  studied.  The  theore- 
tical results  are  compared  with  interferometric  measurements  obtained  in  the 
UTIAS  Hypervelocity  Shock  Tube  for  argon  boundary  layers  on  a flat  plate  and 
on  the  shock- tube  sidewall  behind  a shock  wave  under  close  initial  conditions. 

In  Chapter  2,  the  basic  equations  for  laminar  boundary-layer  flows 
in  partially-ionized  monatomic  gases  are  discussed  and  transformed.  The  basic 
assumptions  are  evaluated.  The  transport  properties  and  chemical-reaction  rates 
are  considered  using  the  known  elastic  and  inelastic- scattering  cross-sections 
for  an  argon  plasma  ( see  Chapter  3)  The  initial  and  boundary  conditions  are 
given  in  Chapter  k.  The  implicit  six-point  method  of  the  finite- difference 
scheme  is  presented  and  discussed  in  Chapter  5*  The  analytical  and  experimental 
results  are  compared  in  Chapters  6 and  7 for  flat-plate  and  shock- tube  sidewall 
boundary-layer  flows,  respectively.  Discussions  and  conclusions  are  given  in 
Chapter  8.  The  explanation  of  the  computer  program  is  presented  in  Appendix  A. 

The  computer  program  is  listed  in  Appendix  P . 

1.3  Basic  Assumptions 

In  the  present  analysis,  the  following  basic  assumptions  axe  used. 

(l)  For  a mixture  of  atoms,  ions  and  electrons,  it  will  be  assumed  that  each 
species  has  a Maxwellian-velocity  distribution  with  an  appropriate  temperature. 
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The  governing  equations  for  the  plasma  motion  can  be  obtained  from 
the  Boltzmann  equation  by  employing  approximations  to  the  distribution  functions. 
One  of  the  cases  for  which  the  equations  are  solvable  is  when  the  particles  have 
a Maxwellian- velocity  distribution.  The  assumption  of  a Maxwellian- velocity 
distribution  can  be  justified  when  all  the  gradients  in  the  macroscopic  properties 
of  the  plasma  are  small  and  no  external  forces  act  on  the  plasma . In  such  a case 
the  plasma  flow  is  an  isentropic  flow,  and  it  can  be  proven  that  the  velocity 
distribution  is  Maxwellian.  The  condition  that  the  velocity  distribution  for 
the  electrons  and  ions  is  near  a Maxwellian  distribution  is  that  the  Larmor 
radius  is  much  greater  than  the  mean-free-path,  or  the  elastic- collision  frequency 
is  to  be  large  during  the  time- evolution  process. 

Under  this  assumption  the  evaluation  of  the  binary-collision  integrals 
in  the  macroscopic  equations  can  be  greatly  simplified.  This  assumption  should 
be  valid  for  the  region  of  the  boundary  layer  except  the  sheath  region  adjacent 
to  the  wall  where  the  electron  and  ion- number  densities  are  very  low. 

(2)  Only  a singly  ionized  species  is  considered.  The  electron-number  density 
can  be  assumed  equal  to  the  iori-number  density.  The  plasma  is  quasi-neutral. 
Therefore,  the  effects  of  elastic  and  magnetic  fields  on  the  boundary- layer 
structure  are  neglected.  The  essential  requirement  for  quasi-charge  neutrality 
is  that  the  Debye  length  is  much  smaller  than  the  characteristic  length  of  the 
problem  (Ref.  18) . The  ambipolar  character  of  the  diffusion  process  results 
from  this  assumption  providing  that  no  electric  currents  cross  the  boundary. 

The  temperatures  of  heavy  particles  and  electrons  considered  here  are 
much  smaller  than  the  ionization  temperature.  Therefore  the  assumption  that 
ions  are  singly  ionized  is  valid  in  general.  However,  in  the  region  adjacent 
to  the  wall,  a space-charge  sheath  exists  wherein  the  gas  is  no  longer  quasi- 
neutral. The  she’ath  is  composed  of  excess  ions,  yielding  an  electric  field  to 
repel  electrons.  Therefore  special  consideration  of  the  sheath  region  is  needed. 

(3)  The  atom  and  ion  temperatures  are  equal.  Therefore,  atoms  and  ions  have 
the  same  velocity.  This  assumption  can.be  justified  since  the  mass  of  the  ions 
is  almost  equal,  to  that  of  atoms,  and  therefore  only  few  collisions  between  atoms 
and  ions  axe  necessary  to  reach  a common  temperature.  Jaffrin  (Ref.  27)  has 
investigated  the  structure  of  a steady  plane  shock  in  a partially  ionizing  gas 
using  the  Navier-Stokes  equations.  He  showed  that  ion  temperature  is  almost 
identical  with  atom  temperature  in  the  whole  region.  However,  the  collisional 
energy-transfer  processes  between  electrons  and  heavy  particles  are  relatively 
slow,  giving  rise  to  a situation  that  electrons  may  have  a temperature  much 
different  from  that  of  heavy  particles.  It  is  shown  from  the  analysis  of  shock 
structure  that  electrons  have  a much  lower  temperature  than  the  heavy  particles 
in  the  ionization- relaxation  zone.  However,  in  a boundary  layer,  electrons  may 
have  a higher  temperature  than  the  heavy  particles  near  the  cool  surfaces. 

Additional  assumptions  made  in  the  present  analysis  are  described  in 
Chapters  2 and  3* 

1.4  Regions  of  Flow  Near  the  Wall 

As  suggested  by  Dix  (Ref.  18),  three  distinct  regions  exist  near  the 
wall:  (l)  Continuum- flow  region;  away  from  the  wall,  the  gas  is  quasi -neutral, 
the  ion-diffusion  velocity  is  small,  and  the  continuum  equations  are  valid.  The 
boundary-layer  equations  described  in  this  report  should  be  applicable.  (2)  Tran- 
sition region;  near  the  wall  but  not  adjacent  to  it,  the  gas  remains  quasi-neutral, 
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but  the  ion-diffusion  velocity  is  of  the  same  order  as  the  ion  sound  velocity. 
The  continuum  equations  are  not  valid.  (3)  Sheath  region;  in  the  region 
adjacent  to  the  wall,  quasi -charge  neutrality  collapses  near  a catalytic 
boundary.  In  this  so-called  sheath  region  with  thickness  on  the  order  of  a 
Debye  length,  strong  electric- field  strengths  are  created  in  order  to  prevent 
the  plasma  from  being  broken  down  in  a very  short  time.  The  sheath  is  composed 
of  excess  ions,  yielding  an  electric  field  to  repel  electrons.  The  Boltzmann 
equation  is  required  in  the  latter  two  regions. 

For  the  flow  conditions  considered  in  the  present  analysis,  the  thick- 
nesses of 'the  transition  and  sheath  regions  are  very  small  compared  with  the 
boundary-layer  thickness.  Since  a major  difficulty  exists  in  the  solution  of 
the  transition  region,  it  is  neglected  and  the  solutions  at  the  edge  of  the 
sheath  region  are  regarded  as  the  wall  conditions  of  the  boundary- layer  flow. 

The  physical  phenomena  of  the  sheath  region  can  be  described  as 
follows.  Whenever  a charged  particle  strikes  an  absorbing  surface,  this 
particle  loses  its  charge  by  recombination  on  the  surface.  Thus,  solid 
surfaces  act  as  sinks  for  charged  particles.  Electrons  have  much  larger 
thermal  velocities  than  the  ions.  Consequently,  per  unit  time  more  electrons 
are  likely  to  strike  the  surface  than  the  slower  ions.  As  the  electrons 
diffuse  in  the  general  direction  of  the  surface,  the  slow  ions  retard  the 
diffusion  by  setting  up  an  electrostatic  field.  This  process  is  called  amibi- 
polar  diffusion,  and  the  associated  electric-potential  field  falls  in  the 
direction  of  the  charge  motion.  Immediately  next  to  the  wall,  the  electron- 
number  density  becomes  too  low  to  carry  the  ions,  and  the  potential  of  the 
surface  and  the  ion-diffusion  motion  take  over.  Therefore,  a sheath  of  high 
electric  field  exists.  Two  methods  can  be  applied  to  the  analysis  of  the 
sheath  region:  (l)  the  continuum- sheath  theory,  and  (2)  the  colli sionless- 
sheath  theory.  In  the  present  analysis,  the  simpler  method  of  colli sionless- 
sheath  theory  is  considered  and  described  in  Section  4.1. 


2.  BASIC  EQUATIONS  AMD  TRANSFORMATION 

2.1  Boundary  Layer  Equations  for  Ionizing  Monatomic  Gases 

A partially  ionized  monatomic  gas  or  plasma  is  considered  consisting 
of  a mixture  of  atoms,  ions  and  electrons.  For  each  species  the  macroscopic 
balance  equations  can  be  expressed  by  using  the  plasma  macroscopic  properties, 
as  shown  below  (Ref.  19), 

§t  tns^  £3  + tns  ^s  ^>3  = I(^s)  (1) 

where  the  quantity  <0S>  is  the  average  of  the  property  0a,  ns  is  the  number 
density  of  species  s,  l(0s)  is  the  source  term  of  property  0S,  and  Vs  is  the 
total  velocity  of  a particle  of  species  s.  The  source  term  expresses  the  change 
in  <0g>  as  a result  of  both  external  influence  (i.e.,  electric,  magnetic  and 
gravitational  fields)  and  internal  influence  (i.e.,  chemical  reactions,  heat 
transfer,  radiation,  diffusion  and  viscosity) . The  continuity,  momentum  and 
energy  equations  are  obtained  by  putting  0S  = ms,  msVs^  and  l/2  msVs«3vs<3  + €int 
respectively,  where  m is  the  mass  of  particle  and  eint  is  the  internal  energy 
of  particle. 
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The  thermodynamic  quantities  used  in  this  work  are  presented  in 
Appendix  B. 

The  general  formulation  of  the  gas dynamic  conservation  equations 
for  individual  species  in  a nonequilibrium  partially  ionized  gas  mixture  has 
been  discussed  by  Appleton  and  Bray  (Ref.  20),  Kaufman  (Ref.  21).  Grewal  and 
Talbot  (Ref.  22),  Spitzer  (Ref.  23)  and  Igra  (Ref.  19) . 

Following  the  above  considerations  the  basic  equations  for  a boundary- 
layer  flow  of  a partially  ionized  gas  are  given  by 


Continuity  equation  for  plasma: 

Sc  (pU)  + |y  (**)  = 0 


Momentum  equation  for  plasma: 

Energy  equation  for  plasma: 


pu  i + pv  i = - b (qc  + V + h (“ u § ) ■ Sr 


(2) 


(3) 


Conservation  equation  for  electron  species: 


pu  S + pT  § = - 1?  + v< 


(5) 


Energy  equation  for  electrons: 

S^e  ^e  ^e  dPe  * 

sr  ♦ sr  = " sr  ♦ * 5T  - i 


+ q. 
ce  Me 


) + V 


cjpe 

i ST  + Qel  + W 
(6) 


with  u,  v as  the  velocities  of  the  plasma  in  x,  y directions,  x coordinate 
along  the  body  surface  and  y normal  to  it;  p,  plasma  density;  p,  plasma  pressure: 
£»  Plas»a-viscosity  coefficient;  H,  total  enthalpy  of  the  plasma;  q~,  plasma 
heat- conduction  flux;  qa,  plasma  diffusive- energy  flux;  Qr,  plasma  radiation- 
energy  loss;  a,  degree  of  ionization;  Vi,  ion  diffusion  velocity  relative  to 
v;  ma,  mass  of  atom  (or  ion);  ne,  electron-nuntoer  density  production  rate;  h-, 
electron  specific  enthalpy  defined  in  Bq.  (7);  Qel>  rate  of  thermal  energy 
given  to  free  electrons  by  elastic  collisions;  Qinei,  inelastic  energy- transfer 
rate;  pe,  electron  pressure;  subscript  e denotes  electron  encounter.  The 
following  quantities  are  applied  in  Eqs.  (2) -(6): 


p = ma(na  + ne) 
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» = (na  + neVa  + "e^e 
H = C (Ta  + OTe)  + RT-j-Ot  + u2/2 


a = V(na  + ne) 

ar  ar 
qc  = • + V sr  - \ sr 


*d  = p.V^RT,  H-  Cl) 


PiVi  = " **>e 


da  . a(l  - a)  a f Te  M 
1“  S(Te/Ta)  ^ ^ 


h = ac  T 

e P e 


CP  -1R 


k - y»a 


^ce  ^ 


&T. 


*de  = CPWi 
pe  = nekBTe 


(7) 
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where  subscripts  a,  i and  e denote  atom,  ion  and  electron,  respectively;  T, 
temperature;  Tj,  ionization  temperature;  n,  number  density;  R,  gas  constant; 

X,  thermal  conductivity  coefficient,  Da,  anbipolar-diffusion  coefficient;  kg, 
Boltzmann  constant . 

The  rate  of  radiant  energy  loss  Qr  of  a plasma  consists  of  the  rates 
of  energy  loss  by  continuum  radiation  and  by  line  radiation.  In  order  to 
simplify  the  calculation,  the  rate  of  the  line  radiant-energy  loss  for  the 
argon  plasma  is  assumed  equal  to  its  continuum  radiation-energy  loss.  This 
assumption  has  been  discussed  in  Refs.  24  and  25.  Based  on  the  assumption  of 
local- temperature  equilibrium,  Qr  is  given  by 


= e6  ne 

3*76  m^/2  c3  h*^  */C^ 


(hvc  + Ve>«  Zeff 


(8) 


with  e,  electron  charge;  h,  Planck  constant;  c,  speed  of  light;  vc,  cut-off 
frequency;  g,  Gaunt  factor;  Zeff,  effective  nuclear  charge. 
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It  has  been  shown  that  excitation  to  the  first  state  is  rate- 
controlling for  the  overall  ionization  process.  We  assume  that  atoms  in  the 
ground  level  are  excited  to  the  first-excitation  level  by  collision  with 
other  particles,  then  excited  atoms  are  ionized  by  subsequent  collisions. 

The  rates  of  reaction  among  levels  higher  than  the  first  are  assumed  to  be 
in  thermal  equilibrium  with  the  electrons  in  the  entire  flow . The  following 
reactions  are  considered  for  the  collisional-ianization  processes: 


fa  + 

Ar  + Ar  2 Ar  + Ar  + e 


ra 


(9) 


fe 

Ar  + e Ar  + e + e 


re 


with  kf  and  kr  as  forward  and  backward-rate  coefficients. 

Based  on  the  two- tender ature  two-step  model  of  Hoffert  and  Lien 
(Ref . 26) , the  electron- nunber  density  production  rate  ne  can  be  written  as 


ne  = (Ae)a  + (ne)e 


(10) 


where  (ng)a  and  (ne)e  denote  the  net  electron  nunber  density  production  rates 
by  atom-atom  collision  and  electron-atom  collision,  respectively.  The  following 
equations  are  used  for  the  electron  nunber  density  production  rates: 


• 2 2 
(n)  =k~n  -k  nn 
' e'a  fa  a ea  a e 

(n)  = k„  n n - k n ^ 

' e'e  fe  a e re  e 


(10a) 


The  elastic  energy- transfer  rate  Qei  is  the  sum  of  the  rates  of 
thermal  energy  given  to  the  free  electrons  by  electron- atom  and  electron-ion 
elastic  collisions  (Ref.  27): 


V = 3n< 


(l) 


(vea  + ,’ei>VTa  ‘ Te> 


(11) 


where  vea  and  vai  are  elastic- collision  frequencies  due  to  electron-atom  and 
electron-ion  encounters,  respectively. 

The  inelastic  energy- transfer  rate  Qinel  is  the  sum  of  the  rates  of 
thermal  energy  given  to  the  free-electrans  by  electron-atom  and  electron-ion- 
electron  inelastic  collisions  and  by  bremsstrahlung . The  latter  is  neglected 
in  the  boundary- layer  flow  since  it  is  small  compared  with  the  former.  For 
the  two-step  model,  Qinel  is  given  by, 


r 


^inel  " 


where  the  term  for  the  creation  energy  of  electron  due  to  atom- atom  ionization 
collisions  is  very  small  and  can  be  neglected  in  inviscid  and  viscous  flows. 

In  order  to  simplify  the  present  analysis , two  approximations  are 
made  for  the  boundary- layer  flows: 

(1)  PlVt  - - pDa  g 


(2)  v • Vpe  ~ 0,  V±  ^ ~ 0 

Approximation  (l)  has  been  widely  accepted  by  many  authors  (for  example, 

Refs.  12  and  28  in  the  analysis  of  two- temper attire  boundary-layer  flows  in 
ionizing  gases.  Approximation  (2)  has  been  used  by  Chung  and  Mullen  (Ref. 

28)  and  Takano  and  Akamatsu  (Ref.  l4)  since  these  terms  are  very  small 
compared  with  others  on  the  RHS  of  Eq.  (6).  However,  in  the  analysis  of 
inviscid  flow  for  ionizing  gases  (for  example,  the  analysis  of  shock-wave 
structure  given  in  Refs.  24  and  25,  approximation  (2)  should  not  be  made. 

Under  these  approximations,  the  total  energy  equation,  Eq.  (4),  for 
the  plasma  becomes 

<’uSE  + 0'r5-35(* 

+ h f + CpTe>  | ] + h (““!)-%  <l4> 

where  X = Xa  + Xi  • 

The  conservation  equation  for  electron  species,  Eq.  (5),  becomes 


pU  & + ^ SF  “ Sr  V.  ^a  ^ ) + man« 


written  as 


Using  Eq.  (15),  the  electron- energy  equation,  Eq.  (6),  can  be  re- 


cpa  pusf  + 


P',5T  ] = If  ( \ 3jr)  * ^ ^ 'W 


- V - <Vi + 1 V.><V.  - 1 V.<V. 
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The  total- enthalpy  equation,  Eq.  (l4) , can  be  rewritten  in  terms  of 
Ta  by  using  Eqs.  (3),  (15)  and  (16): 


0p[pu^5+fW5r]  = ,,i  + |F(^^s)+t‘(^) 

- S>1  - VA).  ■ «R  (17) 

The  basic  equations  for  the  boundary- layer  flows  are  given  by  Eqs. 
(2),  (3),  (15),  (l6)  and  (17)  with  five  unknowns:  u,  v,  a,  Ta  and  Te.  The 
boundary  conditions  of  these  equations  are  discussed  in  Chapter  4. 

2-2  Transformation  of  Boundary-Layer  Equations 

The  similarity  transformation  coordinates  are  applied: 


Sr, 


X 

5(x)  '/  Wb*1 

o 

u r 

Tj(x,y)  = -2-  pdy 
kpT  J 


V 

where  the  subscript  6 denotes  the  edge  of  the  boundary  layer. 
From  Eq.  (l8), 

H - 

to  “ P^6U6 

aa_pu6 


(18) 


(19) 


By  enploying  the  transformed  continuity  equation, 

and  the  transformed  convective  operator, 

ou  3^  = dj  [ (•**  at  a \ ,3  1 

puS;+p',3y  2S  a*  I 25V  Si  3T  3s  35  J ‘ fSi  | 

the  basic  equations  are  transformed  from  the  coordinates  (x,y)  to  the  coor- 
dinates ({,q).  Here, 
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f\[ 

o 6 

The  transformed  equations  for  momentum,  electron  species,  atom 
temperature  and  electron  temperature  are 


[Cf"]’  + ff"  + 0f 


r fs  - f2 


■ 2S  [ f’  If-  -It  f"  ] <20) 


[fc-I 


+ fz'  - p zf  + 
z 


We 


, 2 2 

fs0'  ] +f0'  + Frr0f"2  - ci-7'' 

l J p a6  a r>  ao  K 


2 

i 

p'-aS 


p5^6b 


» 

pT«'  1 + fw®'  + Vf®'  - ®I  Vf'8 

e J e 


+ ^el  " (*BTI  + 2 " 2 WV 


W» 


pC  T _ 
^p  e6 


’ 2«V  [ f'  H - If  ®'  ] <23> 

where  the  prime  denotes  d/cft)  and  the  following  definitions  are  used: 


f' 


z 

e 

e 


u 

C «-BL 

U5 

P5P5 

- 2L 

Sc  = -£r- 

a 

6 

<*a 

T 

fiC 

_ a 

Pr  = — £ 

Ta6 

T 

[iC 

e 

Pr  = t-2- 

T = 
e6 

e \ 

(24) 

Contd. . 
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6 = ^ 
pz  ag  d| 


2|  ^aS 
Ta6  d* 


(24) 

Contd. 


We  note  that  solutions  of  Eqs.  (20)- (23)  are  strongly  dependent  on 
the  thermal  properties  of  the  ionizing  gas.  The  transport  propc.  ies  are 
calculated  from  the  elastic- scattering  cross-sections  and  the  chemical- 
reaction-rate  coefficients  are  calculated  from  the  inelastic- scattering 
cross-sections  for  the  ionizing  gas.  In  the  following  chapter,  the  thermal 
properties  of  ionizing  argon  are  discussed. 


3 • THERMAL  PROPERTIES  OF  AN  ARGON  PLASMA. 


3-1  Elastic-Scattering  Cross-Sections 

The  elastic- scattering  cross-sections  are  used  in  determining  the 
transport  properties  of  ionizing  gases.  They  will  be  evaluated  here  from 
experimental  results  for  argon.  The  average  atom-atom  elastic- collision 
cross-section  craa  is  obtained  from  the  values  of  the  viscosity  coefficient 
given  by  Amdur  and  Mason  (Ref . 29)  • At  high  temperature  n » 31  x 10"7  t&3/4 
g/cm-sec,  which  corresponds  to 


oaa  = 1.7  x lO"14/^’25  cm2  (25) 

Experimental  data  complied  by  Fay  (Ref.  30)  show  that  the  average 
atom-ion  elastic  cross-section  cr^  is  much  bigger  than  the  atom-atom  elastic 
cross-section  because  of  the  charge-exchange  mechanism.  This  cross-section 
decreases  very  slowly  with  the  atom  (or  ion)  temperature  and  will  be  taken  as 

or^  = 2.454  x lO-14/^’09  cm2  (26) 


The  average  momentum- transfer  cross-section  between  electrons  and 
atoms  aea  for  argon  was  calculated  by  Devoto  (Ref.  3l)  using  the  momentum- 
transfer  cross-section  determined  by  Frost  and  Phelps  (Ref.  32).  An  approxi- 
mate value  of  (rea  for  argon  by  curve  fitting  is 

(0.713  - 4.5  x 10"4T  + 1.5  x 10"7T  2)  x 10"16  cm2 

e e 

for  Te  < 3000  K (27) 

(-0.488  + 3-96  x IQ"4!)  x 10_l6  cm2  for  T > 3000  K 

The  average  momentum- transfer  cross-section  between  electrons  and 
electrons  (ree  can  be  obtained  by  assuming  the  relative  kinetic  energy  of 
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electrons  equal  to  1.5  kgTe  in  the  Coulomb  scattering  cross-section: 


Note  that  if  the  above  assumption  does  not  apply,  one  may  obtain  the  following 
form  by  using  the  Maxwellian  distribution  in  the  aeneenad  Coulomb-  scattering 
cross-section: 


Similarly,  the  average  elastic- scattering  cross-section  between 
ions  and  ions,  Oii  is  given  by 


Since  Te/me  » Ta/ma,  the  electron  temperature  is  the  relevant 
temperature  in  the  calculation  of  ion-electron  c ollision  cross-section, 
therefore , 


ei 


ee 


(31) 


3.2  Transport  Properties  of  Ionizing  Argon 

The  kinetic  theory  of  gases  provides  a means  of  estimating  the 
transport  coefficients  of  a partially-ionized  gas.  In  this  section, 
transport  properties  of  partially-ionized  argon  are  considered  as  based  on 
the  mixture  rule  of  Fay  and  Kenp  (Ref.  5). 


The  viscosity  of  plasma  can  be  calculated  as 


1 + 
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ai 


1 + 
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where 


u 


/ Va 
V W“a 


(32) 


is  the  mean  thermal  speed  of  the  atoms.  The  electrons  make  no  contribution 
to  the  viscosity  because  of  their  extremely  low  mass  compared  with  atoms 
and  ions. 


The  ambipolar  diffusion  coefficient  Da  is  defined  in  terms  of  the 
atom-ion  diffusion  coefficient  Dai  by  • 


1 + a ai 


(33) 


or  approximately  as 


( cm2/sec) 


m 


where  the  contribution  of  the  electron  temperature  due  to  electron-ion 
collision  is  negligible  owing  to  the  small  electron  mass. 


The  thermal  conduction  coefficients  for  atoms,  ions  and  electrons 
may  be  written  as  (Ref.  27) 
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(35) 

(36) 

(37) 


3 -3  Inelastic-Scattering  Cross-Sections  and  Reaction-Rate  Coefficients 

I 

The  forward- rate  coefficient  kf,  Eq.(lO),  can  be  obtained  from 
kinetic  theory  by  computing  the  collision  rate  between  two  particles . The 
calculation  requires  knowledge  of  the  dependence  of  the  inelastic- colli si on 
cross-sections  for  the  first  and  higher  excitation  steps  on  incident  energy. 
For  a two-step  model  considered  here,  a knowledge  of  the  energy  dependence  of 
the  cross-section  for  the  first  excited  state  is  required.  Since  only  the 
energy  dependence  of  inelastic  cross-sections  near  the  threshold  energy  is 
jjnportant  in  the  calculation  of  the  rate  coefficients,  same  knowledge  about 
the  energy  dependence  of  inelastic  cross-sections  near  the  threshold  is 
required.  The  inelastic  cross-sections  thus  obtained  are  often  expressed  in 
terms  of  the  following  relation: 


**(*)  * °o  (X  - f ) 

where  Oq  is  a constant  and  e*  is  the  threshold  energy. 
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Nevertheless,  by  making  use  of  Wigner's  R-matrix  theory,  Eu  and 
Liu  (Ref.  33)  have  obtained  a general  form  for  leading  threshold  behaviour 
of  inelastic  cross-sections  in  the  form 

°*(e)  = °-0  (x  - f)7 

which  fits  the  experimental  data  fairly  well  within  the  experimental  errors. 

For  the  present  analysis,  a reasonably  good  approximation  to  this 
cross-section  is  given  by  the  linear  relationship: 


oJt»(  e)  = S^e  - e*)  for  e > e* 


(38) 


where  e is  the  kinetic  energy  (in  centre- of -mass  coordinates)  of  particle  b 
(b  can  be  atcm  or  electron),  e*  is  the  first  excitation  energy  of  particle  a, 
and  Sat  is  first- excitation  collision  cross-section  slope. 

The  forward  rate  coefficient  appearing  in  Eq.  (10)  can  be  written  as, 


wv  ■ s*ab  [ f ( (i^ + x)  ) 

where  T*  is  the  first  excitation  temperature  of  particle  a,  and  this  rate 
coefficient  must  be  divided  by  two  for  like- like  particle  collisions. 


From  a comparison  of  theoretical  and  experimental  results  for  argon 
shock-wave  structure,  we  found  (Ref.  24)  that  Saa  = 1.0  x 10-19  cm2/eV.  A 
more  recent  electron-atom  excitation  cross-section  constant  Sne  * 4.9  x 10” 18 
cn^/eV  for  argon  obtained  by  Zapesochnyi  and  Felston  (Ref.  34)  is  used  here. 
Therefore,  kfa  and  kfe  yield: 

kfa(Ta)  = 1*4  x 10-2V*5(^  + 2 ^ e (cm3/sec)  (4o) 

, \ -T*  /T 

kfe(Te)  = 2-63  x 10-l6Tg-5^  + 2 J e e (cm3/ sec)  (4l) 

Hoffert  and  Lien  (Ref.  26)  used  a chemical  equilibrium  concept  for 
the  present  chemical-nonequilibrium  case  to  determine  kra  and  kre.  However, 
for  low  temperatures  these  results  are  not  valid  and  the  ionic-recombination 
theory  based  on  the  classical  electron- inpact  cross-section  is  needed.  In 
order  to  avoid  the  difficulty  of  determination  of  the  reverse  reaction-rate 
coefficients,  a critical  temperature  Tc  is  defined  which  separates  the  high 
and  low- temperature  regions.  This  critical  tenperature  can  be  obtained  by 
ensuring  the  continuity  of  the  rate  coefficients  at  Tc.  For  the  electron- 
catalyzed  reactions,  Hinnov  and  Hirschberg  (Ref.  35)  have  obtained  an  empirical 
relation  for  the  reverse  reaction-rate  coefficient  at  low  temperature  (Ta  < 
4000  K) . The  following  reverse  reaction-rate  coefficient  kre  for  electron- 
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The  physical  meaning  of  the  cut-off  of  kra  at  low  temperature  is 
that  at  low  Ta  the  reverse  atom-atom  reaction  rate  is  frozen  at  some  particular 
rate  and  the  re-excitation  from  the  first  excited  state  is  not  rate- controlling 
for  atom-atom  collisions.  In  general,  the  reaction  rates  due  to  atom-atom 
collisions  are  very  small  compared  with  those  due  to  electron-atom  collisions 
for  atom  temperature  below  about  15,000  K.  Therefore,  the  atom- catalyzed 
reactions  can  be  neglected  for  Ta  < 15,000  K,  in  a flat-plate  boundary- layer 
analysis  where  the  flow  has  cooled  significantly.  However,  for  the  case  of 
a shock-tube  sidewall  boundary  layer  near  the  shock  front  where  the  atom 
temperature  is  large  (~  25,000  K),  atom- catalyzed  reactions  are  more  important 
than  the  electron- catalyzed  reactions  and  kra  must  be  retained.  Byron  et  al 
(Ref.  36)  have  shown  that  for  the  low-temperature  case,  de-excitation  from 
other  than  the  first-excited  state  can  be  rate  controlling  in  the  recombination 
process.  For  the  present  two-step  model,  the  approximation  made  in  Eq.  (43)  is 
necessary  in  order  to  avoid  the  unknown  physical  effects  due  to  a very  large 
value  of  kra*  It  is  also  worth  noting  that  a large  value  of  kra  destabilizes 
the  finite-difference  scheme. 

Another  method  for  evaluating  the  rate  of  atom- catalyzed  reactions 
is  to  cut  off  kra  and  to  limit  the  rate  of  recombination  reaction  of  atom-atom 
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collision  at  some  particular  value.  The  following  equation  can  be  applied: 


(n)  = k„  n 2 - k n n 2 if  (n  ) > 0 

v e'a  fa  a ra  a e ^ e'a 

= 0 ^ (Ae)a  < 0 

4.  BOUNDARY  AMD  INITIAL  CONDITIONS  FOR  BOUNDARY-LAYER  FLOWS 
4.1  Boundary  Conditions 

The  boundary  conditions  for  Eqs.  (2)-(6)  are 


(44) 


y = 0: 


u =uw 
v = 0 

T = T (or  ST  /Sy  = 0 for  zero  heat  transfer) 
8*  w a 


u 

U6 

a 

11 
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Cn 

T 

= T 

a6 

a 

T 

= T e 

e 

e6 

(45) 


where  the  wall  values  uw  and  Tw  are  usually  given.  The  other  values  ug,  ag, 
Tag  and  Teg  are  determined  from  the  inviscid-flow  region  (see  Section  4.3). 


The  boundary  conditions,  Eq.  (45),  Tor  the  transformed  equations, 
Eqs.  (20) -(23),  are 

ij  = 0:  f = 0 

f'  - Vu6 

0 = Tw/Tag  (or  0'  = 0 for  zero  heat  transfer) 

(46) 

T)  -»  00:  f ' = 1 

Z = 1 


0=1 
e = 1 
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The  other  required  boundary  conditions  are  z and  8 at  the  wall; 
these  can  be  obtained  frcm  the  following  nethods : 

(l)  Wall-Sheath  Properties : 

For  the  case  in  which  the  wall  is  at  floating  potential,  the  first 
equation  is  obtained  from  Langrauir-probe  theory  for  a Maxwellian  distribution 
of  electrons: 


n <V  > 
es  e 


( ‘ ^ ) ' n«eTi " 0 


where  <Vg>  = (Q^^/m*)1/2 , V±  = (k^ Jv^)1/2,  £4>  is  the  potential  difference 
between  wall  and  plasma  and  subscript  s denotes  the  value  evaluated  at  the 
sheath  edge. 

The  second  equation  is  the  continuity  of  mass  flow  of  ions  at  the 
outer  edge  of  the  sheath: 


psDas  (10  = PsasVi 


The  energy  equation  at  the  edge  of  the  sheath  is 


T ^e  1 

*B  W " [p0!Ws  = (2kBTes  + 

-IS 

where  Va  is  the  ambipolar-diffusion  velocity. 
From  Eqs . (47)  and  (48)  we  obtain 


n <7  > 
es  e 


(■#) »’ 


and  Eq.  (49)  becomes 


«•  . (so)  sS  life  01/2 

TJ  ' ' it  n 


W PwU6  WW 


where  Vi6  = (l^Te6/ma)1/2  ^ - in(ma/21jme)1/2ew. 

This  model,  based  on  continuity  at  the  sheath  edge,  was  widely  used 
by  many  authors,  for  example,  Camac  and  Kemp  (Ref.  37),  Dix  (Ref.  18),  Nishida 
and  Natsuoka  (Ref.  12),  Sherman  et  al  (Ref.  ll)  and  Mansfeld  (Ref.  15).  However, 
Mansfeld  (Ref.  15)  mentioned  that  the  artificial  boundary  condition  used  for 
the  two-temperature  equilibrium  model  leads  to  values  of  ne  near  the  wall  which 
seem  to  be  in  much  better  agreement  with  the  experimental  results  than  the 
values  obtained  from  an  electric- sheath  consideration.  He  concluded  that  the 
validity  of  the  boundary  condition  for  z and  0 at  the  wall  derived  from  a 
presently  incomplete  description  of  the  sheath  is  still  unknown. 
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(2)  Catalytic  Wall  Model 


By  analogy  with  the  dissociated  boundary- layer  flow,  the  wall  is 
assumed  to  be  a catalytic  wall  when  an  equilibrium  composition  is  used  at 
the  wall.  For  a one-tenqperature  equilibrium  model  where  Taw  « 300  K,  the 
boundary  condition  for  z at  the  wall  is  approximately  given  by 


z » 0 
w 


For  a two- teiqper at ure  equilibrium  catalytic  wall  model,  z*  is  given  by 


z 


w 


= g 


where  g is  a constant  value. 

For  a cooled-wall  case  with  a wall  sheath  model,  lakano  and 
Akamatsu  (Ref.  l4)  have  shown  that 


zw  ~ 0 (10’2)A££ 
e;  ~ 0 (10"1*) 

where  Re  is  the  Reynolds  number. 

We  also  note  that  Nishida  and  Matsuoka  (Ref.  12)  have  shewn  that 
the  slope  of  the  electron  temperature  at  the  wall  is  almost  equal  to  aero. 
Man sf eld  (Ref.  15)  has  obtained  the  following  results  even  when  X is  vary 
small: 


zw  « 0 (52a) 

% ~ 0 (52b) 


In  order  to  obtain  better  agreement  between  theory  and  experiment, 
Tseng  and  Talbot  (Ref.  13)  have  used  a measured  value  of  Zw  as  the  wall 
boundary  condition: 


z = 0 .02 
w 


4.2  Conqpatlbility  Conditions 

At  the  edge  of  the  boundary  layer,  the  following  boundary  conditions 
must  be  satisfied: 
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f"  = 0 


z'  =0 


0'  - 0 


0'  = 0 


By  using  the  above  boundary  conditions,  the  compatibility  conditions  at 
i|  -+eo  can  be  obtained  from  Eqs.  (20)-(23): 


p = -J\  ■-  JL1 

Ws  Pb<x* 
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These  conditions  must  be  satisfied  in  the  calculations  in  order 
to  avoid  a discontinuity  in  the  gradients  of  the  dependent  variables  at  the 
edge  of  the  boundary  layer. 

The  value  dpg/dt  appearing  in  Eq.  (5^)  must  be  obtained  either 
from  experiment  or  theory.  The  following  considerations  should  be  noted 
in  the  calculation: 

(1)  If  the  external  effects  (for  example,  an  unsteady  effect)  or 
interactions  (for  example , the  interactions  with  a shock  wave  or  an  expansion 
wave)  occur  in  the  inviscid-flow  region,  the  value  dpg/d£  is  obtained  from 
the  solution  of  the  inviscid  flow  with  these  effects  or  interactions  taken 
into  account. 

(2)  If  there  is  no  external  effect  or  interaction  in  the  inviscid 
flow,  another  equation  is  needed  to  form  a complete  set  of  equations  with 
five  unknowns:  dpg/dj,  0f,  pz,  Pra  and  Pre.  This  equation  is  obtained  from 
the  continuity  equation  of  the  plasma  flow,  Eq.  (l) , 
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After  same  algebraic  arrangement, 

Pf  = A#u  2 ^Ta6PT&  + a6Te6PTe  + a6Te6Pz^ 
where  0ja,  Pre  and  |3Z  are  calculated  from  Eq.  (5*0  and 


It  was  shown  by  Blottner  (Ref.  17)  that  a swallowing  of  the 
inviscid  flow  into  the  boundary  layer  is  necessary  in  order  to  satisfy  the 
compatibility  conditions. 

Since  same  approximations  have  been  made  in  the  boundary-layer 
equations  with  respect  to  the  inviscid- flow  equations,  the  values  dp 5/di, 

Pf > Pz>  0To  and  PTe  obtained  from  the  above  equations  should  be  slightly 
different  from  that  obtained  from  the  inviscid- flow  region.  If  there  is  no 
approximation  made  in  the  boundary-layer  equations,  or  the  inviscid- flow 
equations  are  obtained  directly  by  letting  all  d/ciy  terms  equal  to  zero  in 
the  boundary-layer  equations,  then  the  above  method  provides  the  same 
results  that  would  be  obtained  from  the  inviscid- flow  region. 

A local- similarity  method  was  applied  by  Brown  and  Mitchner  (Ref. 
38)  in  predicting  the  electron- temperature  and  electron-number- density 
profiles  of  a flat-plate  boundary-layer  plasma.  They  predicted  that  the 
electron  temperature  at  the  edge  of  the  boundary  layer  is  smaller  than  the 
atom  temperature  and  explained  it  as  due  to  the  radiation-energy  loss. 
However,  it  is  clear  that  the  compatibility  conditions  described  above 
were  not  applied  in  their  calculations.  The  electron  temperature  at  the 
edge  of  the  boundary  layer  must  be  calculated  from  the  equations  for  the 
inviscid  flow  and  not  from  the  boundary-layer  equations.  In  their  calcula- 
tion, in  order  to  satisfy  the  boundary  conditions  at  the  edge  of  the  boundary 
layer,  Eq.  (53)j  the  values  of  the  degree  of  ionization  and  the  electron 
temperature  at  the  edge  were  adjusted. 
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The  solutions  for  Ug,  Tag,  Teg,  ag  and.  pg  must  be  obtained  from 
the  equations  for  the  inviscid-flow  region.  The  quasi-one-dimensional 
equations  for  inviscid  flow  are  obtained  from  Eqs.  (l)-(5)  by  letting 

d/dy  = 0: 


Ie  (0“>  - 0 

(55) 

pu  5“  + |E  = o 

M dx  dx 

(56) 

dH 

‘"‘Sc  "-«R 

(57) 

da 

pu  — = m n 
dx  a e 

(58) 

$E  <°Ve>  ‘ “ ST  * «» l + V 

(59) 

Equations  (55) -(59)  have  been  solved  by  Glass  and  Liu  (Ref.  24) 
for  the  shock-wave  structure  of  ionizing  argon  and  by  Glass  et  al  (Ref.  25) 
for  a krypton  shock-wave  structure.  The  inviscid  flow  generated  by  a 
shock  wave  can  be  separated  into  two  zones:  (l)  an  ionization- relaxation 
zone  and  (2)  a radiative- cooling  zone.  In  the  relaxation  zone  the  elastic 
and  inelastic-collision  processes  are  important  while  in  the  radiative- 
cooling  zone  the  radiation- energy  loss  is  significant.  Equations  (55) -(59) 
provide  a unified  treatment  applicable  to  both  zones.  However,  from  our 
numerical  experience  in  solving  the  shock-wave  structure,  a complete  solution 

for  the  radiative- cooling  region  requires  a small  step-size  to  be  stable.  As 
the  plasma  is  nearly  in  equilibrium,  values  for  ug,  Tag,  Teg,  ag  and  pg  in 
the  cooling  region  can  also  be  obtained  approximately  by  solving  only  Eqs . 
(55)-(57)  together  with  the  Saha  equation.  Whitten  (Ref.  39)  has  shown 
that  the  error  in  using  a radiant  equilibrium  model  is  within  2$  of  the 
present  nonequilibrium  model. 

4 .4  Initial  Conditions 

The  initial  profiles  are  required  for  a finite- difference  method. 

At  the  start  of  the  boundary  layer,  | = 0,  and  the  partial-differential 
equations  become  ordinary-differential  equations.  At  £ = 0:  0f  = 0ra  = 

0Te  = Pz  = Oj  the  following  ordinary-differential  equations  are  obtained 
from  Eqs.  (20)-(23): 


[Cf"]*  + ff"  =0  (60) 
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Equations  (6o)-(63)  with  two-point  boundary  conditions  can  be 
solved  using  the  usual  iteration  techniques  (for  example,  the  Newton-Raphscn 
method).  A subroutine  BEGIN  for  solving  Eqs.  (6o)-(63)  is  discussed  and 
presented  in  Appendix  F . 


5.  FINIIE-D1FFKHENCE  METHOD 


5.1  Mathematical  Considerations 


Numerical  methods  have  developed  rapidly  in  the  last  decade  and 
solutions  have  now  been  found  to  many  systems  of  simultaneous  equations 
which,  prior  to  the  development  of  the  digital  computer,  could  not  be  solved 
because  of  the  immense  amount  of  calculation  required.  One  of  the  most 
common  numerical  techniques  for  solving  partial-differential  equations  is 
the  finite-difference  method,  where  the  differential  equations  are  replaced 
by  a large  number  of  difference  equations,  which  are  then  solved  by  various 
algebraic  methods. 

The  means  of  solving  simultaneous  algebraic  equations  can  be 
divided  into  two  general  types  called  direct  and  indirect  methods.  Direct 
methods,  which  include  elimination  and  matrix  inversion  techniques,  require 
a finite  number  of  steps  to  obtain  an  exact  solution.  Indirect  methods 
theoretically  require  an  infinite  number  of  steps  to  obtain  a solution  but 
often  can  provide  a sufficiently  accurate  solution  in  a much  smaller  number 
of  steps  than  would  be  required  with  a direct  method.  The  large  number  of 
difference  equations  resulting  from  partial-differential  equations  make 
direct  methods  impractical  for  solving  a problem.  For  indirect  methods, 
for  example,  the  modified-Leibmann  method,  or  over-relaxation  method,  are 
free  from  round-off  errors  and  have  the  additional  advantage  that  they  can 
often  be  adapted  to  solve  nonlinear  equations.  The  direct  methods  are 
usually  applied  in  the  parabolic  or  hyperbolic- type  partial  differential 
equations  having  sides  with  an  open  boundary.  The  indirect  methods  are 
applied  to  elliptic- type  equations  with  a closed  boundary. 

The  finite- difference  method  for  linear  partial-differential 
equations  has  been  well  established.  Unfortunately,  methods  for  solving 
nonlinear  algebraic  equations  are  lagging  far  behind.  Recently,  because 
of  the  large  nunber  of  physical  and  engineering  problems  which  are  described 
by  nonlinear  equations  and  the  prospects  which  the  computer  offers  for  their 
solutions,  the  techniques  for  solving  nonlinear  algebraic  equations  has 
become  an  active  field  of  mathematical  research. 

Two  general  methods  have  been  developed  for  the  solution  of  a set 
of  simultaneous  nonlinear  algebraic  equations.  The  first  is  called  by 
Greenspan  (Ref.  ho)  the  nonlinear-Liebmann  method,  which  involves  linear- 
izing the  equations  by  putting  known  values  into  the  nonlinear  terms  and 
requires  iteration.  The  resultant  set  of  linear  algebraic  equations  is  then 
solved  by  the  extrapolated-Liebmann  method  or  some  other  method.  The  process 
of  iteration  is  continued  until  all  the  residuals  are  suitably  small.  The 
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second  method,  known  either  as  the  generalized  Newton's  method  or  nonlinear 
over- relaxation , is  an  iterative  procedure  where  the  grid  is  scanned  in 
order.  The  generalized  Newton's  method  is  the  faster  of  the  two  nonlinear 
methods  according  to  Greenspan,  hut  this  is  based  only  on  examples  and  not 
a general  mathematical  theory. 

The  above  iteration  techniques  are  general  mathematical  methods 
developed  to  be  applicable  to  a wide  range  of  nonlinear  differential  equations. 
When  the  generalized  Newton's  method  was  applied  to  same  engineering  problems, 
it  was  found  to  work  satisfactorily  only  for  one-dimensional  cases. 

When  a nonlinear  differential  equation  is  a description  of  some 
physical  situation,  the  nature  of  the  nonlinearity  is  known  and  often  a 
special  numerical  method  can  be  devised  to  control  the  nonlinearity  during 
a relaxation- type  iteration  procedure.  When  the  generalized  Newton's  method 
has  serious  drawbacks  then  the  above  becomes  necessary.  For  example,  the 
projection  method  (Ref.  4l)  can  be  applied  in  order  to  control  the  nonlinear 
terms. 

For  the  nonlinear  parabolic  differential  equations,  the  projection 
method  is  usually  applied  in  order  to  control  the  nonlinear  term.  For  example, 
we  consider  the  following  nonlinear  equation: 

/ x S^F  _ ^F 

a(F)  7T  ‘ -St 
dx 

The  finite  difference  analogs  used  in  solving  this  simple  type  of  equation 
are  centred  around  the  time  level  t + &t/2,  and  the  coefficient  a(F)  must 
be  evaluated  at  this  time  level.  The  simplest  method  of  solving  this  type 
of  equation  is  an  iteration  process  where  for  all  grid  points  i,  the  terms 
a(Fi)  are  first  evaluated  using  the  values  of  Fi  at  time  t.  Substitution 
of  these  terms  into  the  difference  equations  results  in  a set  of  linear 
equations  which  are  easily  solved  for  the  function  F at  the  time  t + 6t. 

The  coefficients  a(Fj_)  are  then  re-evaluated,  using  for  F^  the  average  of 
its  value  at  time  t and  the  newly  calculated  value  for  time  t + &t.  After 
substitution  of  the  newly  calculated  a(Fj.)  terms,  the  difference  equations 
are  again  solved  for  the  function  Fi  at  time  t + 5t.  This  iteration  is 
repeated  until  the  function  Fi  determined  in  two  successive  iterations  agree 
within  a predetermined  tolerance.  The  nonlinear  terms  have  been  projected 
forward  to  the  level  t + 6t/2.  More  sophisticated  techniques  use  a Taylor 
series  in  conjunction  with  the  finite-difference  analog  of  the  original 
equation  to  project  the  nonlinear  terms  to  the  half-time  level.  The  above 
projection  method  can  be  used  in  the  usual  boundary-layer  equations.  Further 
discussion  on  the  nonlinear  partial- differential  equations  in  engineering 
applications  can  be  found  in  Ref.  42. 

Many  methods  can  be  applied  for  the  parabolic- type  partial- 
differential  equations,  for  example,  the  explicit  method,  implicit  method, 
Crank-Nicolson  method,  DuFort-Frankel  method,  Saul'yev  method,  and  the 
explicit  and  implicit  alternating- direction  methods.  In  the  explicit  method, 
usually  undesirable  restrictions  on  the  step- size  increment  occur  in  the 
computation.  The  implicit  method  can  overcome  this  difficulty  at  the  expense 
of  a somewhat  more  complicated  calculation  procedure.  However,  the  discreti- 
zation errors  for  both  methods  are  still  too  large.  The  Crank-Nicolson 
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method  is  of  the  implicit  type  and  also  can  decrease  the  discretization 
errors.  It  is  always  stable  and  the  error  is  of  second  order. 


The  solution  of  equations  resulting  from  the  implicit  or  Crank- 
Nicolson  method  can  be  obtained  by  any  elimination  technique.  However, 
since  the  resulting  equations  have  the  form  of  the  tridiagonal-type  the 
complete  algorithm  method,  which  has  less  of  an  iteration  scheme  than  the 
Gaussian  or  Gauss-Seidal  elimination  methods,  can  be  applied.  A general 
form  which  can  be  applied  by  explicit,  implicit  or  Crank- Nicolson  method  is 
developed  in  the  present  work.  An  excellent  review  for  the  finite  differ- 
ence method  of  solution  of  the  boundary-layer  equations  has  been  given  by 
Blottner  (Ref.  17) . 

5 -2  Finite  Difference  Equations 

The  nonlinear  equations,  Eqs.  (20)-(23),  with  the  boundary 
conditions  of  a mixed  Neumann/Dirichlet  type  are  solved  numerically  by 
the  finite-difference  method.  An  implicit  six-point  finite-difference 
scheme  is  applied. 

These  equations  are  first  linearized  in  a form  suitable  for  an 
iteration  scheme.  Blottner  (Ref.  17)  has  stated  that  the  order  of  the 
equations  is  important.  The  momentum  equation  is  solved  first,  and  the 
species  must  be  solved  before  the  atom  temperature.  Therefore,  these 
linearized  equations  can  be  written  in  order  as  follows: 
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The  superscript  p denotes  the  order  of  the  iteration  process  and 
the  quantities  without  the  superscript  denote  those  evaluated  at  the  p-1 
iteration  order.  F = df/<h)  or 


f = J Fdi,  (68) 

o 

These  linearized  equations  are  of  the  second  order  and  are  solved 
for  the  unknowns  F,  z,  0 and  9 in  that  order.  The  derivatives  and  the 
integral  in  the  rj-direction  are  then  expressed  by  three-point  difference 
formulae.  The  derivatives  in  the  |-direction  are  approximated  by  a forward- 
difference  scheme . 

Let  i and  j be  the  indices  of  the  tj,  |- coordinates  for  the 
difference  net  at  the  point  considered  in  Fig.  1.  Any  function  W (F,  z,  6 
or  9)  is  written  in  terms  of  the  values  of  two  adjacent  points  in  the 
{-direction  as 


W = AW(i,  j+i)  + (i  . X)w(i,  j) 

where  X is  a weighting  factor  which  can  be  suitably  adjusted  for  improving 
the  convergence  of  the  iteration  scheme: 


0;  Explicit  method 


Crank-Nicolson  method 


1:  Implicit  method 


In  this  formulation  either  equal  intervals  or  nonequal  intervals 
in  ij-direction  can  be  used.  In  the  present  case,  the  interval  in  ^-direction 
is  increased  in  a geometric  progression  as 
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where  k is  a constant  which  is  set  with  a value  slightly  greater  than  unity. 
The  following  derivatives  are  applied: 
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where  Arj^  is  the  first  interval  in  rj-direction. 

The  quantity  g is  evaluated  at  a point  between  two  adjacent 

points  as 

5 = * *d+3^1_^*d 


The  values  f and  fg  at  a (i,j+l)  point  are  given  by 
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where 


?5U,  d+D  - W>  d+D  - *(*,  J)3 

The  linearized  equations,  Eqs.  (64) -(67),  nay  be  written  in  the 
following  common  form  with 

+ x^wj^  + x^v/1*  = x^V^  + x^  (74) 

where  = P,  = z,  W^  = e and  = 0.  The  expressions  Xj^  (i  = 1 
to  4)  for  the  momentum,  species,  atcm  temperature  and  electron  temperature 
equations  are  listed  in  Appendix  C . 

Substituting  Eqs.  (69) -(71)  into  Eq.  (74),  the  following  equation 
is  obtained: 

A^fc)W(t)(i-l,  d+1)  +B^tVt)(i,  J+l)  + c[fcVfc)(i+l,  j+1)  = (75) 

where  t = 1 to  4;  i = 1 to  N and  j = 1 to  M. 
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The  boundary  conditions  for  the  finite  difference  equations  are: 


w(1)(o,  j+1)  = ~ 

6 

v/ 2)(0 , j+1)  =0 

for  z =0 
w 

= W(2)(l,  j+1) 

• Vi 

for  z;  = ^ 

W<3)(0>  j+1)  = ~- 

a6 

W(4)(0,  j+1)  = W(4)(l,  j+1) 

- 

for  e;  = 4, 

= W(4)(l,  j+1) 

for  0'  = 0 
w 

i = N (l)  -»eo) 


W^^(N,  j+1)  = 


W(2)(N,  j+1) 


= 1 


W^(N,  j+1) 


= 1 


W^(N,  j+1)  = 


The  maximum  value  of  i (or  N)  which  represents  the  freestream 
condition  can  be  determined  as  follows:  After  the  values  k and  are 
chosen  ( see  next  section) , the  results  for  all  grid  points  at  a fixed  j 
with  an  arbitrary  value  N are  compared  with  that  calculated  from  the  N+l 
value.  If  they  do  not  have  the  same  value,  then  N is  increased  until  the 
results  of  using  N and  N+l  have  identical  values.  In  order  to  guarantee 
that  the  value  N used  represents  the  freestream  condition  for  all  j , N is  given 
by  N = No  + 20,  where  N0  is  the  minimum  value  of  N at  x = 0.  The  maYinnrm  value 
of  tj  is  obtained  from  the  following  equation, 

N 

w 1 

i=l 


The  computational  scheme  is  an  iterative  one.  The  f 
momentum  equation  is  first  solved  with  assumed  distributions  of  species  and 
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atom  and  electron  temperatures.  The  resulting  velocity  field  is  employed  for 
the  species  equation.  The  resulting  species  field  is  then  applied  to  the 
equations  for  atom  and  electron  temperatures,  respectively.  The  new  species 
and  atom  and  electron  temperature  distributions  are  then  used  to  replace  the 
assumed  one  and  the  process  is  continued  until  the  solutions  converge  to 
satisfy  a preset  criterion. 

The  convergence  criterion  of  the  system  of  difference  equations  to 
the  differential  equations  has  not  been  investigated.  However,  Douglas  et  al 
(Ref.  43)  stated  that  an  implicit  or  Crank-Nicolson  difference  scheme  is 
convergent  for  an  equation  of  the  type  given  by  Eq.  (74) . 

In  order  to  avoid  third-order  derivatives  in  the  momentum  equation, 
Blottner  (Ref.  7)  has  introduced  a transformed  normal  velocity  and  retained 
the  continuity  equation.  However,  in  the  present  method,  the  stream  function 
is  introduced  and  the  momentum  equation  is  written  as  a second-order  equation. 
The  partial  differential  equation  for  momentum  involves  f and  fg.  The  value 
of  f can  be  readily  obtained  from  an  integration  once  the  value  of  the  tangen- 
tial velocity  component  across  the  boundary  layer  is  known.  The  same  method 
was  applied  by  Fannelop  (Ref.  44).  Sells  (Ref.  45)  also  used  the  same  implicit 
finite-difference  method  for  a laminar  conpressible  boundary  layer  and  Chan 
(Ref.  46)  for  a turbulent  incompressible  boundary  layer. 

5 *3  Accuracy  and  Stability 


The  accuracy  of  the  numerical  solution  has  to  be  better  than  the 
accuracy  of  the  different  physical  models.  The  models  used  for  the  descrip- 
tion of  transport  properties,  chemical  reactions  and  sheath  theory  are  not 
supposed  to  have  a higher  accuracy  than  0(10"^).  The  accuracy  of  the  experi- 
mental results  is  at  best  0(10_2) . Therefore,  it  seems  sufficient  to  achieve 
numerical  results  which  are  accurate  to  0(10“^). 


The  accuracy  of  the  numerical  method  can  be  achieved  to  any  small 
order,  say  0(l0-5) , at  the  expense  of  computation  time.  Once  the  accuracy  of 
the  problem  is  determined,  the  upper  bounds  for  A|  and  A»|  are  posed. 

The  accuracy  of  any  numerical  method  can  be  checked  by:  (l)  varying 
the  | and  tj  increments,  (2)  disturbing  the  boundary  conditions  slightly,  (3) 
applying  the  difference-differential  method  or  other  numerical  methods,  (4) 
applying  a stability  analysis  to  the  linearized  equations,  and  (5)  applying 
the  method  to  a simple  problem  that  can  be  solved  analytically.  In  the 
present  analysis  items  (l)  to  (4)  were  applied  to  check  the  accuracy  of  the 
Crank-Nicolson  scheme.  The  step  sizes  used  in  the  calculation  are  decreased 
by  half  and  the  solutions  follow  this  in  a way  corresponding  to  the  orders 
of  the  local  truncation  errors  and  remain  stable.  A small  disturbance  to  the 
input  data  gives  a small  change  in  the  solution.  The  results  obtained  by 
using  X = 0.5  ( Crank-Nicolson  method)  and  X = 1 (inplicit  method)  were 
compared  and  the  accuracy  was  within  0(l0“2) . 

In  the  present  analysis,  nonequidistant  step  sizes  were  used  in 
order  to  decrease  the  computation  time.  However,  the  nonequidistant  discreti- 
zation may  lead  to  larger  inaccuracy  and  can  even  spoil  the  solution  completely. 
Also  the  determination  of  error  bounds  for  nonequidistant  step  sizes  is  more 
complicated  than  an  equidistant  regular  network.  Fortunately,  the  nonequidistant 
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discretization  can  be  checked  in  the  present  program  by  coinparing  the  results 
with  k = 1 and  with  k / 1.  It  was  shown  that  the  nonequidistant  discretization 
does  not  lead  to  larger  inaccuracy  or  spoil  the  solution  in  the  present  program. 

As  mentioned  in  Section  5.2,  the  stability  criteria  for  a system 
of  nonlinear  partial-differential  equations  are  difficult  to  determine.  When 
the  Crank-Nicolson  scheme  is  applied,  bounded  oscillations  in  5-direction 
appear  in  the  analysis.  Crandall  (Ref.  47)  showed  that  at  a relatively  large 
step  size  in  |-direction,  bounded  oscillations  are  possible  even  for  linear 
equations.  According  to  Smolderen  (Ref.  48)  this  may  be  even  worse  in  the 
case  of  nonlinear  equations.  Douglas  (Ref.  49)  also  mentioned  the  possibility 
of  oscillations  in  the  solution  using  a Crank-Nicolson  scheme  with  boundary 
conditions  of  the  mixed  type.  These  oscillations  do  not  occur  when  a backward- 
implicit  method  is  used.  However,  a smaller  step  size,  A5 , then  the  one  used 
in  the  Crank-Nicolson  scheme  is  needed  for  the  backward- implicit  scheme  for 
the  same  accuracy.  This  results  in  more  computation  time.  For  the  present 
case,  the  oscillation  can  be  controlled  by  a suitable  choice  of  the  parameter 
X.  From  our  experience  of  the  present  analysis,  the  Crank-Nicolson  scheme 
provides  a bound  oscillation  in  5-direction.  This  oscillation  does  not  riamp 
out  even  for  a very  small  step-size  A|.  Therefore,  decreasing  step  size  A| 
is  not  the  best  way  to  get  more  accuracy.  This  oscillation  can  be  checked 
by  using  different  values  of  X in  the  calculations  starting  from  an  optimum 
value  0.5  to  a maximum  value  of  1.0,  and  the  smallest  value  of  A in  this 
range  which  just  eliminates  the  oscillations  is  the  one  to  use.  In  the 
present  analysis,  the  best  value  of  X was  found  to  be  0.75  where  the  oscilla- 
tion damps  out  in  the  first  few  steps.  With  X = 0.75j  the  discretization 
error  is  expected  to  be  0(A5^"5),  which  is  smaller  than  0(A|)  for  a backward- 
implicit  scheme. 

For  exanple  (see  Chapter  6,  for  the  flat-plate  boundary  layer) 

U5  - 3-53  x 10?  cm/ sec,  Ta6  = Te6  = 1.059  x 101*  K,  a6  = 0.031  and  ps  = 1200 
torr  was  used  for  all  x.  The  step- size  A|  was  increased  with  x.  At  first 
the  value  1( of  thetweight  parameter  X = 0.5  was  used,  and  in  all  runs  oscilla- 
tion in  fw  and  started  at  the  first  step  downstream.  This  oscillation 
tended  to  be  small  at  i * 45.  In  another  run  X = O.75  was  used,  and  this 
time  the  oscillation  damped  out  within  two  steps  downstream  and  did  not  re- 
appear. Figure  2 shows  how  the  values  fw  and  zw  oscillate  as  the  step  number 
j increases.  Consequently,  the  Crank-Nicolson  scheme  (X  = 0.5)  was  abandoned 
and  X - 0.75  was  used  in  the  present  analysis. 

Attenpts  were  made  to  relate  this  oscillation  to  errors  due  to  the 
step-size  Atj  and  A|  with  the  aid  of  formulae  like  the  Richardson- extrapolation 
rule,  but  a satisfactory  answer  was  not  obtained.  Therefore,  the  oscillations 
were  not  induced  only  by  a finite  step-size  as  they  were  considerably  larger 
than  the  errors  expected  from  such  step  sizes. 

5.4  Transformation  of  Coordinates 

The  conditions  at  the  edge  of  boundary  layer,  Ug,  <25,  Ta§,  Teg  and 
P&j  resulting  from  the  solutions  of  the  inviscid  flow  equations  (see  Section 
4.3)  are  a function  of  x.  A table  of  these  edge  properties  as  a function  of 
x was  used  and  the  interpolation  was  applied  to  obtain  the  edge  conditions 
for  any  value  of  x.  However,  the  finite  difference  procedure  was  applied  for 
the  transformed  ( |,tj) -coordinates . Therefore,  the  final  results  must  be 
related  back  to  the  physical  (x,y) -coordinates. 


The  relation  between  the  transformed  coordinate  | and  the  distance 
x along  the  surface  can  be  determined  from  Eq.  (18) . The  relation  between 
£ and  x can  be  determined  by  the  following  equation: 
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where  the  three-point  formula  for  nonequidistant  step  sizes  applied  in  Eq. 
(72)  can  be  used  for  Eq.  (76). 

6.  FLAT-PLATE  BOUNDARY -LAYER  FLOWS  IN  IONIZING  ARGON 
6.1  General  Considerations 


Experiments  on  shock-wave  structure  and  the  boundary-layer  flows 
induced  by  a strong  shock  wave  were  recently  conducted  at  UTIAS  in  the 
Hypervelocity  Shock  Tube.  These  experiments  provided  unique  and  reliable 
interferometric  data  for  both  types  of  flows.  Experiments  on  shock-wave 
structure  were  conducted  by  Bristow  (Ref.  50) , Brimelow  (Ref.  51),  Tang 
(Ref.  52)  and  Whitten  (Ref.  39 )>  and  on  the  flat-plate  and  the  shock- tube 
sidewall  boundary  layers  by  Whitten  (Ref.  39)  and  Brimelow  (Ref.  51), 
respectively.  Comparisons  of  numerical  and  experimental  results  on  shock- 
wave  structures  are  given  in  Refs.  24  and  25  for  argon  and  krypton,  respec- 
tively. 


Measurements  of  ionizing  flat-plate  boundary- layer  flows  have  been 
reported  by  same  authors  (Tseng  and  Talbot,  Ref.  13;  Brown  and  Mitchener, 

Ref.  38;  Bredfeldt  et  al.  Ref.  53)  for  low  temperatures  and  low  electron- 
nuntoer  densities.  Under  these  conditions,  the  radiation- energy  loss  in  the 
plasma  is  small  and  can  be  neglected  for  both  the  inviscid  and  viscous  flow 
regions . Thus,  the  conditions  at  the  edge  of  the  boundary  layer  can  generally 
be  calculated  frcm  a nonradiant  model  in  which  the  free  stream- flow  quantities 
are  constant  and  independent  of  time. 

Figure  3 shows  schematically  the  experimental  generation  of  a flat- 
plate  boundary-layer  flow  over  sin  airfoil  model  with  a sharp  expansion  corner 
in  the  UTIAS  10  cm  x 18  cm  Hypervelocity  Shock  Tube.  Such  a boundary  layer 
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can  be  regarded  as  developing  in  a steady  flow  if  the  shock  wave  is  travelling 
at  a constant  velocity  and  the  radiant  energy  loss  is  small.  However,  for  the 
case  of  a flow  induced  by  a stronger  shock  wave  (Me  > 13) , the  radiation- 
energy  loss  becomes  significant  and  the  boundary  layer  develops  in  a somewhat 
unsteady  (nonurJJbrm)  flow. 

A typical  analysis  of  the  inviscid  radiant  flow  behind  a shock  wave 
travelling  at  constant  velocity  is  shown  in  Fig.  4.  In  this  case,  the  shock 
wave  is  moving  at  a Mach  number  Ms  = 12.8  into  quiescent  argon  at  pQ  = 5.01 
torr  and  To  = 297  K and  has  been  shown  at  a location  xs  = 40  cm  past  the 
leading  edge  of  the  flat  plate.-  The  gradients  result  from  the  radiation- 
energy  loss.  Clearly,  as  the  nonstationary  shock  wave  travels  along  the  tube, 
the  inviscid- flow  conditions  above  the  plate  change  with  time,  introducing  an 
unsteady  (nonuniform)  effect.  Along  the  flat  plate,  it  is  seen  that  ug  and  p6 
decrease  as  x increases  from  the  leading  edge  (or  as  one  moves  closer  to  the 
shock  front)  while  Tag,*  Teg  and  ccg  all  increase.  However,  the  fact  that  the 
inviscid  flow  with  respect  to  the  plate  is  unsteady  rather  than  steady  with 
flow  gradients,  is  emphasized  by  a consideration  of  the  overall  momentum 
equation  for  the  plasma.  In  a steady,  one- dimensional  inviscid  flow, 

dug  dpg 

P6U5  dx  " dx 

That  is,  the  velocity  gradient  and  pressure  gradient  should  have  opposite 
signs,  which  is  not  the  case  in  this  flow. 

It  can  be  seen,  however,  that  the  variations  of  ug,  Tag*  Teg  and  pg 
are  quite  small,  particularly  as  the  distance,  Xs,  from  the  leading  edge  to 
the  shock  wave  increases.  Under  these  circumstances,  it  is  reasonable  to 
regard  the  flat-plate  boundary-layer  flow  as  quasi- steady  to  a first-order 
approximation , such  that  the  steady-flow  analysis  described  previously  can  be 
applied.  It  should  be  mentioned  that  the  relative  changes  in  ag  are  slightly 
larger,  however,  and  the  full  extent  of  this  effect  on  the  quasi- steady- flow 
assumption  is  not  known  at  the  present  time. 

For  comparisons  of  boundary- layer  profiles  measured  at  a position 
Xffl  with  analytical  predictions,  the  inviscid  flow  conditions  at  xm  were 
assumed  to  prevail  over  the  entire  freestream  region  (all  x)  in  order  to 
satisfy  the  steady- flow  assumption  in  the  analysis.  The  initial  conditions 
for  the  shock  wave  and  the  freestream  quantities  resulting  from  the  radiant 
inviscid  flow  are  listed  in  Table  1. 

In  the  finite  difference  analysis,  the  best  value  of  the  weight 
parameter  was  found  to  be  A = 0.75,  where  the  oscillations  damped  out  during 
the  first  few  steps.  Case  2 of  Table  1 was  run  with  k = 1.05,  Ajj  = 0.05  and 
N = 46  by  using  a step-size  Ax  started  with  0.01  cm  at  x = 0 and  increased 
to  0.2  cm  at  x = l4  cm.  At  first  the  value^of  the  weight  parameter  A = 0.5 
was  used,  and  in  all  runs  oscillations  in  fw  and  started  at  the  first-step 
downstream.  The  oscillations  were  bounded  and  tended  to  become  smftii  as 
x increased.  The  other  run  with  A = 0.75  was  used,  and  this  time  the  oscilla- 
tions damped  out  within  two  downstream  steps  and  did  not  reappear.  Therefore 
the  Crank- Nicolson  scheme  (A  = 0.5)  was  abandoned  and  A = 0.75  was  used  in  the 
analysis. 





6.2  Comparison  of  Theoretical  and  Experimental  Results 

For  the  case  of  Mg  = l6.6  and.  pQ  = 4.8l  torr,  the  initial  flow 
profiles  (at  x = 0)  are  shown  in  Fig.  5*  The  nondimensional  electron- 
temperature  profile  8 is  unity  for  all  tj,  resulting  from  the  electric- sheath 
condition  8'  = 0 at  x = 0.  The  atom- temperature  profile  9 increases  from 
0.029  at  tj  = 0 to  1 at  tj  = 1.95,  and  reaches  a maximum  value  1.04  at  tj  = 2.4 
and  then  approaches  unity  at  tj  = 4.25.  The  normalized  velocity  profile  f' 
increases  from  zero  at  tj  = 0 to  1 at  tj  = 4.5,  while  the  normalized  degree  of 
ionization  profile  z increases  from  0 at  tj  = 0 to  1 at  tj  = 3*5* 

The  variations  of  the  transport  properties  at  x = 0 of  Pr,  Sc,  C 
and  Pre  with  tj  are  shown  in  Fig.  6.  The  Prandtl  number  for  the  heavy 
particles  Pr  is  constant  (Pr  = 0.667)  from  its  definition  in  this  analysis. 
The  ratio  of  density-viscosity  product  C decreases  from  2.8  at  tj  = 0 to  1 
at  tj  = 2.23.  The  Schmidt  nuntoer  Sc  increases  from  1.5  at  tj  = 0 to  2.42  at 
Tj  = 1.4,  and  then  decreases  to  the  freestream  value  2.39  at  tj  = 3*12. 
Similarly,  Pre  increases  from  0.035  at  tj  = 0 to  0.52  at  tj  = 2.4  and  then 
decreases  to  the  freestream  value  0.507  at  tj  = 4.79.  These  transport- 
property  parameters  are  functions  of  na,  ne,  Ta,  Te  and  P5.  The  variations 
of  these  parameters  with  tj  have  some  effects  on  the  boundary-layer  structure. 
The  effect  of  Pre  on  the  electron-tenperature  profile  is  more  significant 
than  that  of  C on  the  velocity  profile  and  of  Pr  on  the  atom-temperature 
profile.  The  total-Prandtl  nuntoer  Pr  of  the  plasma  can  be  obtained  from  the 
following  equation: 


The  variations  of  the  flow  profiles  with  distance  x are  shown  in 
Fig.  7.  The  velocity  profile  f'  is  almost  independent  of  x.  The  variation 
of  atom  temperature  ratio  9 with  x is  also  small.  Therefore,  the  momentum 
and  atom-temperature  equations  can  be  obtained  approximately  from  a similarity 
assumption.  However,  significant  variations  of  the  degree  of  ionization  ratio 
z and  electron-tenperature  ratio  9 with  x do  occur,  as  shown,  and  errors  will 
result  from  a similarity  assumption.  The  degree  of  ionization  a and  the 
electron  temperature  Te  decrease  as  x increases  at  a constant  tj.  At  tj  = 0, 

0w  decreases  from  1 at  x = 0 to  O.87  at  x = l4  cm. 

The  variations  of  Pre  with  tj  for  x = 0 and  x = l4  cm  are  shown  in 
Fig.  8.  It  is  seen  that  for  x > 0,  Pre  exceeds  the  freestream  value  (0.507), 
up  to  tj  ~ 2.5. 


In  addition  to  the  profiles  of  the  various  flow  quantities  across 
the  boundary  layer,  parameters  that  characterize  the  skin  friction,  heat 
transfer  due  to  conduction  and  diffusion  processes  and  thickness  of  ttye  boundary 
layer  are  important . The  variations  of  the  skin  friction  parameter  fw  and  the 
heat  transfer  parameters  9^  and  for  conduction  and  diffusion  processes, 
respectively,  with  distance  x are  shown  in  Fig.  9*  The  values  of  f-^  are 
almost  independent  of  x,  while  9^  increases  at  small  x and  approaches  a nearly 
constant  value  for  large  x.  The  quantity  z^  decreases  significantly  for  x < 2 
cm  and  approaches  a nearly  constant  value  for  large  x.  The  boundary-layer- 
displacement  thickness  5*  is  plotted  in  Fig.  10  as  a function  of  x.  For  x 
greater  than  4 cm,..  6*  increases  almost  linearly.  The  physical -boundary- layer 
thicknesses  with  a flow-quantity  ratio  of  O.995  for  velocity  5f,  degree  of 
ionization  6Z,  atom  temperature  6ra  and  electron  temperature  6re  are  plotted 
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in  Fig.  11.  It  can  be  seen  that  the  boundary-layer  (x  > 3 can)  thickness  for 
velocity  is  greater  than  the  other  boundary-layer  thicknesses  for  this  case. 

The  thickness  of  the  electron  thermal  layer,  Sre,  is  thinner  than  the 
velocity  thickness  which  differs  from  the  result  found  by  Honma  and  Komuro 
(Ref.  l6)  for  a sidewall  boundary- layer  flow.  They  showed  that  the  thickness 
of  the  electron  thermal  layer  is  almost  ten  times  that  of  the  velocity  or 
a-ccot- temperature  thicknesses. 

The  large  change  in  the  chemical-reaction  rates  with  temperature 
has  an  important  effect  on  the  boundary-layer  structure  for  the  case  of  a 
large  degree  of  ionization.  Figure  12  shows  the  effects  of  chemical-reaction 
rates  on  boundary-layer  structure.  The  results  for  a frozen  flow  (he  = 0) 
are  compared  with  those  for  a nonequilibrium  flow  at  x = l4  cm  and  Ms  = l6.6. 
The  profiles  of  velocity  f ’ , and  heavy  particle  temperature  9,  hardly  differ 
for  both  cases.  However,  the  profiles  of  electron  temperature  0,  degree  of 
ionization  z,  and  electron  Prandtl  number  Pre,  are  significantly  affected  by 
the  chemical  reactions.  For  a given  ij,  the  electron  temperature  e is  lower 
for  a frozen  flow  than  for  a nonequilibrium  flow,  while  the  reverse  is  true 
for  the  degree  of  ionization  z. 

Comparisons  of  analysis  with  experimental  results  are  shown  in  Figs. 
13  and  l4  for  plasma  density  and  electron-number  density,  respectively.  Better 
agreement  is  obtained  for  the  measured  plasma  density  profile  with  the  frozen- 
flow  analysis.  However,  poor  agreement  with  analysis  is  obtained  for  the 
electron  number-density  profile  with  either  solution.  The  experimental  data 
show  a significant  bump  in  the  ne  profile  which  is  not  predicted  by  either  the 
nonequilibrium  or  the  frozen- flow  analysis . A similar  bump  appears  in  the 
ex~-  .rimental  data  for  the  degree  of  ionization  profile  shown  in  Fig.  15,  while 
the. e is  no  bump  in  the  analytical  result.  This  disagreement  has  not  been 
resolved. 


For  the  second  case  with  Mg  = 12.8,  pD  = 5*01  torr  and  T0  = 297  K, 
the  nonequilibrium  and  frozen-flow  profiles  at  x = l4  cm  are  shown  in  Figs. 

16  and  17,  respectively,  for  comparison.  It  is  also  shown  that  significant 
differences  exist  for  degree  of  ionization  and  electron  temperature  profiles, 
as  predicted  for  Mg  = l6.6  case.  The  analytical  and  experimental  results 
for  plasma  density  and  electron  number  density  profiles  sire  compared  in  Figs. 

18  and  19,  respectively,  while  the  corresponding  degree  of  ionization  profile 
is  shown  in  Fig.  20.  Unlike  case  1,  the  experimental  plasma- density  data 
shows  better  agreement  with  the  nonequilibrium  or  the  calculated  equilibrium 
similarity- solution  profiles,  which  are  very  close.  The  experimental  results 
for  ne  lie  between  the  analytical  nonequilibrium  and  frozen- flow  profiles.  The 
two-temperature  frozen-flow  solution  predicts  a larger  bump  than  that  obtained 
from  the  experiment,  while  no  bump  occurs  in  the  nonequilibrium  or  equilibrium 
profiles.  The  experimental  data  for  the  degree  of  ionization  lie  closer  to 
the  calculated  frozen  profile  rather  than  the  nonequilibrium  profile. 

The  local- similarity  solution  based  on  thermal  and  chemical  equili- 
brium are  also  plotted  in  Figs.  18-20  for  case  2.  It  is  seen  that  equilibrium 
would  not  occur  in  such  boundary- layer  flows.  Of  the  three  models  used  in 
the  analyses  for  ne  and  a,  the  equilibrium  profiles  provide  the  worst  agreement 
with  the  experimental  results.  The  ne  profile  is  the  most  sensitive  indicator 
of  the  state  of  the  boundary  layer.  The  fact  that  the  density  data  agree 
with  the  three  models  used  in  the  analyses  shows  that  it  is  not  a sensitive 
parameter.  Undoubtedly,  measurements  of  electron  and  heavy-particle  temperatures 
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are  desirable  as  they  would  be  sensitive  indicators  of  the  state  of  the 
• boundary  layer . 


The  disagreement  between  theory  and  experiment  for  ne  (or  a)  may 
result  from:  (i)  the  boundary-layer  flow  is  assumed  to  be  quasi-steady  while 
in  the  experiments  it  is  actually  nonstationary  with  time  resulting  from 
radiation  losses  (which  were  only  partially  accounted  for) , and  the  effects 
of  the  sidewall  boundary  layers  (which  were  not  taken  into  account  at  all 
and  have  flow  effects  similar  to  radiation) . These  effects  will  be  more 
pronounced  at  higher  shock  Mach  numbers  where  the  radiation-energy  loss  is 
significant.  For  example,  the  agreement  between  theory  and  experiment  for 
Ms  =16.6  is  worse  than  that  for  Ms  = 12.8  for  the  ne  and  a profiles.  In 
order  to  assume  that  the  boundary-layer  flow  is  quasi- steady,  the  variations 
of  the  flow  quantities  at  the  edge  of  the  boundary  layer  with  distance  were 
neglected  in  the  theory.  This  error  might  have  a significant  effect  on  the 
boundary- layer  structure,  (ii)  The  assumptions  made  in  the  basic  equations 
(such  as  PiVi  *»  -pDa  bx/b y and  v • Vpe  5=3  0) . the  uncertainty  of  the  para- 
meters used  to  describe  the  elastic  and  inelastic- energy- transfer  rates  and 
the  model  used  for  the  radiation- energy  loss,  may  all  affect  the  entire 
boundary-layer  structure.  However,  the  comparison  of  analysis  and  experiment 
for  the  Mg  = 12.8  case  is  quite  good  and  lends  support  that  the  assumptions 
made  in  the  basic  equations  are  reasonable.  Furthermore,  from  the  shock-wave- 
structure  analyses  (Glass  and  Liu,  Ref.  24;  Glass  et  al.  Ref.  25)  the  parameters 
used  for  the  elastic  and  inelastic- energy  transfer  and  the  radiation  model  are 
considered  as  accurate  within  the  limitations  of  present-day  collision  theory. 

As  noted  earlier  same  possible  errors  in  the  present  analysis  may 
result  by  neglecting  the  re-absorption  of  the  radiation-energy  loss  in  the 
freestream  and  the  effects  of  the  sidewall  boundary  layers  on  the  freestream 
flow.  An  exact  solution  to  the  set  of  simultaneous  ordinary  differential 
equations  for  the  freestream  flow  including  re-absorption  would  be  difficult 
since  the  re-absorption  coefficient  is  a function  of  the  complete  structure 
of  the  radiation  cooling  zone.  The  question  would  also  arise  whether  the 
shock  tube  has  a finite  or  an  infinite  optical  depth.  The  Rosseland  mean- 
free-path  for  argon  at  the  freestream  conditions  of  Mg  = 16.6,  pQ  = 4.8l 
torr  and  To  = 296  K (case  1),  is  about  97  cm*  Therefore,  the  freestream 
plasma  is  optically  thin  and  the  re- absorption  energy  should  be  small.  For 
the  present  calculations  it  was  necessary  to  consider  the  worst  case  when 
there  is  no  re-absorption.  The  shock- tube  sidewall  boundary  layer  will  have 
same  effects  on  the  freestream  conditions . In  general , it  would  be  desirable 
to  study  the  present  flat-plate  boundary  layer  after  the  effects  of  the 
shock-tube  sidewall  boundary  layers  on  the  freestream  were  determined.  In 
the  present  analysis  the  freestream  conditions  were  obtained  under  the 
assumption  that  the  flow  was  one- dimensional;  the  role  of  the  sidewall 
boundary-layer  growth  on  the  inviscid  flow  was  not  considered.  Mirels  (Ref. 

54)  has  shown  that  the  flow  between  the  shock  wave  and  contact  surface  in  an 
actual  shock  tube  is  non-uniform  due  to  the  wall  boundary  layer.  Recently, 
Enflmoto  (Ref.  55)  has  studied  the  effects  of  the  boundary- layer  growth  in 
shock  tubes  of  various  cross-sections  on  the  shock-wave  ionization- relaxation 
process  in  argon.  He  used  Mirels ' boundary- layer  theory  for  a perfect  gas 
to  get  some  estimates.  He  showed  that  by  considering  the  sidewall  boundary- 
layer  effects  the  ionization  relaxation  time  was  shortened.  The  inclusion 
of  re-absorption  of  radiation  energy  and  sidewall  boundary-layer  effects  will 
increase  the  degree  of  ionization  in  the  freestream.  Neglecting  these  two 
effects  in  the  analysis  might  alter  the  freestream  conditions.  In  order  to 
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study  the  effects  of  the  free stream  conditions  on  the  boundary-layer  structure 
for  the  Mg  = 16.6  case,  where  the  theoretical  and  experimental  boundary- layer 
thicknesses  for  ne  axe  0.8  nm  and  2 mm,  respectively  (x  = l4  cm),  three 
different  freestream  values  (at  t = 20  (is,  40  (is  and  60  (is  after  the  passage 
of  the  shock  wave)  were  used  for  comparison.  Even  then,  the  theoretical 
boundary- layer  thicknesses  for  ne  are  between  0.8  mm  and  1.2  mm  and  still 
differed  from  the  experimental  results.  It  is  shown  that  the  boundary- layer 
thickness  for  ne  increases  as  the  degree  of  ionization  in  the  freestream 
increases. 


It  is  worth  noting  that  the  error  due  to  neglecting  the  photo- 
ionization in  the  analysis  may  not  be  small.  Near  the  wall,  where  the 
ionization  due  to  electron-atom  collisions  is  small  compared  with  that  in 
the  freestream,  photons  resulting  from  stimulated  emission  and  radiation 
processes  may  have  an  opportunity  to  ionize  atoms.  The  populations  of  the 
photon  flux  should  be  known  before  any  calculation  on  the  photo-ionization 
rate  can  be  done. 

The  effects  of  chemical  reactions  on  boundary- layer  structure  and 
the  sidewall  boundary  layer  on  the  shock  structure  are  discussed  in  Appendices 
D and  E,  respectively. 

Finally,  from  a comparison  of  the  theoretical  results  of  Mansfeld 
(Ref.  15)  and  the  experimental  results  of  Kuiper  (Ref.  56)  for  a thermal 
Rayleigh  boundary  layer,  it  was  shown  that  the  frozen  models  lead  to  values 
obtained  from  a nonequilibrium  model.  No  bump  in  ne  occurred  in  either  theory 
or  experiment  for  a thermal  Rayleigh  boundary  layer.  Details  of  the  results 
described  in  this  chapter  are  given  in  Ref.  57. 


7.  SHOCK-TUBE  SIDEWALL  BOUNDARY- LAYER  FLOWS  IN  IONIZING  ARGON 
7 .1  General  Considerations 


A considerable  amount  of  theoretical  and  experimental  work  has 
been  done  on  the  prediction  of  heat  transfer  and  the  variations  of  transport 
and  thermodynamic  properties  through  the  shock-tube  sidewall  (or  Rayleigh) 
boundary  layers.  Kuiper  (Ref.  56)  made  an  interferometric  study  of  the 
shock-tube  endwall  boundary-layer  flow.  Mansfeld  (Ref.  15)  compared  his 
numerical  results  with  Kuiper' s experimental  data.  However,  except  for 
Brimelow's  results,  there  are  no  experimental  data  on  shock- tube  sidewall 
boundary-layer  flows  for  ionizing  argon. 

From  the  comparison  of  analyses  and  experiments  for  the  flat-plate 
and  thermal  Rayleigh  boundary  layer  flows,  it  has  been  shown  that  electrons 
have  a temperature  which  is  very  different  from  the  heavy  particles  owing  to 
the  slow  colli sional  energy-transfer  processes  between  electrons  and  heavy 
particles.  In  the  inner  part  of  the  boundary  layer  (i.e.,  near  the  wall)  an 
equilibrium  assumption  is  not  valid. 

Shock-tube  sidewall  boundary  layers  for  ionizing  argon  flows  were 
analyzed  by  KnOOs  (Ref . 1)  for  flows  with  thermal  and  chemical  equilibrium,  Honma 
and  Komuro  (Ref.  16)  for  frozen  flows  and  by  Takano  and  Akamatsu  (Ref.  l4)  for  flows 
with  thermal  and  chemical  nonequilibrium.  Radiation  energy  losses  in  both  free- 
stream and  boundary-layer  flows  were  neglected  in  these  analyses.  Kh55s  has  shown 
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that  the  equilibrium  solution  for  the  degree  of  ionization  is  in  good  agreement 
with  experiment  only  near  the  edge  of  the  boundary  layer.  Honma  and  Kamuro 
found  that  the  thickness  of  the  electron  thermal  layer  for  frozen  flow  is 
about  ten  times  the  thickness  of  the  viscous  boundary  layer,  which  is  in 
contrast  to  the  results  for  flat-plate  (described  in  Chapter  6)  and  thermal 
Rayleigh  boundary-layer  flows. 

The  present  chapter  deals  with  conparisons  of  measurements  and 
related  analyses  of  the  toted- density  and  electron-number-density  profiles  in 
shock-tube  sidewall  boundary  layer  flows.  The  shock- wave -structure  model  of 
Glass  and  Liu  (Ref.  24) , which  includes  radiation-energy  losses,  was  used  and 
Brimelow's  interferometric  data  were  accurately  re-evaluated  for  this  purpose . 

A comparison  was  also  made  between  theoretical  and  experimental  plasma  tenpera- 
tures.  Satisfactory  agreement  has  been  obtained  for  shock-tube  sidewall  boundary- 
layer  flows  in  ionizing  argon. 

7.2  Coup ari son  of  Theoretical  and  Experimental  Results 

The  main  differences  between  the  flat-plate  and  shock- tube  sidewall 
boundary  layers  are:  (i)  the  velocity  profile  for  the  flat-plate  case  increases 
from  zero  at  the  wall  to  the  free stream  value  at  the  edge  for  the  sidewall 
boundary  layer,  (ii)  Ionizing-nonequilibrium  phenomena  occur  in  the  freestream 
flow  behind  the  shock  front.  Consequently,  variations  of  freestream  conditions 
for  the  sidewall  boundary  layer  are  significantly  larger  than  for  a flat  plate, 
(iii)  The  sidewall  boundary  layer  induces  significant  changes  in  the  freestream, 
especially  in  the  ionizing  shock  structure  and  beyond  where  the  radiation  losses 
are  large . 

Recently,  Brimelow's  data  were  re-evaluated  and  compared  with  the 
radiant  shock- structure  model  of  Glass  and  Liu.  The  data  show  that  (i)  no 
bump  occurs  in  the  electron-number-density  profile  in  the  sidewall  boundary- 
layer,  unlike  the  flat-plate  case,  and  (ii)  significant  two-dimensional  and 
wall  effects  on  the  freestream  flow  exist  in  the  sidewall  case  at  higher  shock 
Mach  numbers . 

Accurate  interferometric  determinations  of  the  excitational  cross- 
section  constants  for  argon  and  krypton  atom-atom. collisions  in  the  relaxation 
zone  were  made  by  Glass  et  al  over  ranges  of  the  initial  shock  Mach  number  and 
pressure  of  13  < Ms  < 18  and  3«1  < Po  < 5*2  torr,  respectively.  This  shock  wave 
structure  model  was  used  to  determine  the  range  of  freestream  conditions. 

Two  cases  were  studied  for  Mg  = 13*1,  Po  = 5*6  torr  and  Ms  = 15 .9> 

Po  = 5*10  torr.  The  initial  conditions  and  the  freestream  conditions  at  the 
measuring  station  xm  axe  given  in  Table  2.  Figure  21  shows  schematically  the 
shock-tube  sidewall  boundary  layer  behind  a shock  front. 

• _ • 

Figure  22  shows  a plot  of  the  freestream  conditions  for  Mg  =13.1 
and  Po  = 5 .16.  The  boundary  layer  profiles  were  measured  at  the  cascade  front 
where  the  electron  number  density  is  a maximum  and  the  variations  of  the  free- 
stream conditions  are  significant.  At  the  cascade  front,  the  radiation  energy 
loss  rate  is  maximum. 

The  dimensionless  nonequilibri urn- flow  profiles  of  velocity  F, 
degree  of  ionization  a,  atom  temperature  9 and  electron  tenperature  0 
are  shown  in  Fig.  23.  It  is  seen  that  the  velocity  profile  F is  signi- 
ficantly different  from  that  of  the  flat-plate  boundary  layer  (Fig. 
l6) . The  thickness  of  the  electron  thermal-layer  is  much  thinner 
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than  that  predicted,  by  Honma  and  Komuro.  They  indicated  a thickness  almost 
an  order  greater  than  the  present  result. 


A comparison  of  analysis  with  experimental  data  is  shown  in  Figs. 

2k  and  25  for  plasma  density  and  electron  number  density,  respectively. 

Much  better  agreement  is  obtained  between  the  measured  results  and  the 
nonequilibrium  solutions  for  the  sidewall  than  for  the  flat-plate  boundary  layer. 
Better  agreement  was  obtained  between  the  measured  plasma  density  and  electron- 
number-density  profiles  and  the  frozen- flow  analysis  for  the  flat-pla^q,  boundary 
layer.  A bump  appeared  in  the  experimental  data  for  the  electron-nunber-density 
profile  in  the  flat-plate  case.  No  bump  appeared  in  the  sidewall  case.  The 
agreement  with  the  frozen  solution  is  rather  poor  for  the  sidewall  case.  The 
corresponding  profile  for  the  degree  of  ionization  is  shown  in  Fig.  26.  Owing 
to  the  shape  of  the  velocity  profile,  6*  is  now  negative,  as  expected.  The 
boundary- layer  displacement  thickness,  B*,  is  plotted  in  Fig  27.  It  is  seen 
that  the  displacement  thickness  increases  almost  linearly  with  x,  when  x is 
greater  than  2 cm.  A comparison  with  Fig.  10  shoWs  that  ' ' 

the  sidewall  boundary  layer  is  about  an  order  thicker  for  an  equivalent  x. 

For  the  second  case  with  Ms  = 15-9  and  pQ  = 5.1  torr,  the  freestream 
conditions  for  pg,  ne5  and  ag  are  shown  in  Fig.  28  together  with  experimental 
results.  Good  agreement  is  found.  The  dimensionless  nonequilibrium  and 
frozen- flow  profiles  at  x = 18  cm  are  plotted  in  Fig.  29.  Except  for  the 
velocity  profile  F,  these  are  significantly  different  between  the  nonequilibrium 
and  the  frozen-flow  solutions  for  degree  of  ionization,  and  the  atom  and  electron- 
temperature  profiles.  For  the  flat-plate  case,  profiles  of  velocity  and  atom 
temperature  hardly  differ  for  both  cases.  Chemical  reactions  have  a significant 
effect  on  the  degree  of  ionization  profile. 

Comparisons  of  analysis  with  experimental  results  are  shown  in  Figs. 

30  and  31  for  plasma  density  and  electron-number-density  profiles,  respectively. 
Better  agreement  is  obtained  for  the  measured  plasma- density  profile  with  the 
nonequilibrium  analysis.  However,  only  fair  agreement  with  analysis  is  obtained 
for  the  electron  number  density  profile.  The  latter  is  overpredicted  by  analysis, 
which  is  in  contrast  with  the  results  for  case  1.  It  is  also  shown  that  the 
electron- number  density  continued  to  increase  with  distance  from  the  wall  and 
did  not  reach  asymptotic  values.  This  phenomenon  was  also  observed  in  the  shock 
wave  structure  experiments  (see  Fig.  13,  Glass  and  Liu,  Ref.  24).  In  the 
relaxation  region  near  the  shock- tube  wall,  the  electron  cascade  front  moves  in 
towards  the  wall  slowly  at  first  and  then  very  rapidly.  The  reasons  for  this 
premature  ionization  close  to  the  wall  in  the  experimental  shock  wave  structure 
are  far  from  clear.  One  possibility  considered  was  that  a gas-surface  inter- 
action occurred  between  the  argon  plasma  and  the  chromium-plated  steel  shock- 
tube  wall.  However,  two  experiments  were  carried  out  to  try  and  eliminate 
this  possibility  by  changing  the  surface  material.  However,  no  changes  were 
observed.  This  phenomenon  is  much  more  in  evidence  for  the  stronger  shocks. 

For  this  reason,  experimental  results  for  electron-nunber  density  for  Ms  = 15-9 
are  lower  than  that  predicted.  For  the  previous  case  with  Ms  = 13«1>  the  flow 
parameters  reached  their  asymptotic  values  at  the  edge  of  the  boundary  layer. 

The  corresponding  profile  for  the  degree  of  ionization  is  shown  in  Fig.  32  and 
the  displacement  thickness,  5*,  is  plotted  in  Fig.  33  as  a function  of  distance  x. 
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The  determination  of  electron  temperature  of  a shock-heated  argon 
plasma  has  received  wide  attention.  Among  the  commonly  used  plasma  diagnostics 
for  electron  temperature  determination  in  a flowing  gas  plasma  are  Langmuir 
probing  and  microwave  transmission.  The  difficulty  in  measuring  two  temperatures 
of  a plasma  generated  by  stronger  shock-induced  boundary  layers  has  been  recog- 
nized. In  the  present  case,  the  experimental  temperature  of  argon  plasma  is 
determined  from  the  measured  plasma- density  and  electron- number-density  profiles 
by  assuming  thermal  equilibrium  in  the  boundary-layer  flow.  Even  this  assumed 
experimental  temperature  cannot  represent  the  actual  atom  or  electron  tempera- 
ture. However,  same  interesting  features  can  be  found  from  a comparison  of 
the  theoretical  and  experimental  plasma-  temper  at  ure  profiles.  The  experimental 
plasma-temperature  profile  together  with  the  calculated  two-texnperature  profiles 
are  plotted  in  Fig.  3^  for  Ms  = 13*1  and  pQ  = 5.16  torr.  It  is  shown  that  the 
experimental  plasma  temperature  is  close  to  the  electron  temperature  profile 
near  the  wall.  In  the  outer  part  of  the  boundary  layer,  agreement  between  the 
theoretical  and  the  indirect  experimental  results  is  excellent.  A comparison 
of  the  temperature  profiles  for  Mg  = 15-7  and  pQ  = 5-1  torr  is  shown  in  Fig. 

35. 


From  the  foregoing  comparison  of  theory  and  experiment  for  both 
flat-plate  and  sidewall  boundary-layer  flows,  the  following  differences  are 
observed:  (l)  The  experimental  data  and  the  nonequilibrium  analysis  show 
that  there  is  no  bump  in  the  iiq  profile  for  the  sidewall  boundary  layer. 
However,  such  a bump  is  observed  in  the  experimental  flat-plate  boundary 
layer  profile.  The  frozen  solution  predicts  a bump  for  both  layers.  (2) 

Better  agreement  is  obtained  between  experiment  and  the  nonequilibrium  analysis 
in  the  sidewall  case,  while  the  frozen  solutions  agree  better  with  experiment 
than  the  nonequilibrium  solution  for  the  flat-plate  case.  (3)  More  significant 
differences  exist  between  the  analytical  nonequilibrium  and  frozen  plasma- 
density  profiles  for  the  sidewall  boundary  layer,  while  this  difference  is 
sraal 1 for  the  flat-plate  case.  (4)  Even  though  the  predicted  displacement 
thickness  for  the  sidewall  boundary  layer  is  an  order  of  magnitude  greater,  the 
predicted  and  actual  density,  electron-number  density  and  degree  of  ionization 
layers  have  correspondingly  similar  values.  The  depsity  thickness  is  usually 
about  half  the  electron-number  density  thickness.  It  shows  that  the  total 
density  is  not  a sensitive  indicator  of  flow  variations. 

By  and  large  all  analytical  profiles  are  consistent  for  both  types 
of  boundary  layers.  The  fact  that  the  experimental  data  is  in  oetter  agree- 
ment for  the  sidewall  boundary- layer  flow  may  result  from:.,  (l)  The  unsteady 
effects  on  the  sidewall  boundary  layer  close  to  the  shock  front  are  smialler 
than  those  for  the  flat-plate  cases  where  the  shock  wave  is  threefold  the 
distance  away.  (2)  There  is  significant  ionizing  nonequilibrium  in  the  free- 
stream  flow  in  the  sidewall  case  but  the  flat  plate  is  almost  in  a state  of 
radiant  equilibrium.  (3)  The  degree  of  ionization  at  the  points  measured 
for  the  sidewall  case  is  two  to  threefold  larger  than  that  for  the  flat-plate 
cases  and  may  favour  the  agreement  with  the  nonequilibrium  analytical  profiles. 

It  is  worth  noting  that  the  step  size  used  for  the  sidewall 
boundary  layer,  is  much  smaller  than  for  the  flat-plate  case  and  results  in 
increased  computation  time. 

The  assumption  pj_Vj  « - pDa(  da/Sy)  made  in  the  analysis  was  removed 
and  replaced  by 


4l 


Plvi 

in  the  computer  program.  It  is  shown  that  the  effects  on  the  electron-number- 
density  and  plasma- density  profiles  are  very  small.  This  was  done  to  remove 
instabilities  at  higher  shock-Mach  number  for  the  sidewall  boundary  layer  by 
using  the  above  equation. 

• 

Finally,  from  a comparison  of  the  theoretical  results  of  Hutten- 
Man^feld  (Ref.  15)  and  the  experimental  results  of  Kuiper  (Ref.  56)  for  a 
thermal  Rayldigh  boundary  layer,  it  was  shown  that  no  bump  in  1%  occurred  in 
both  results.  It  should  be  noted  that  in  Kuiper' s experiment  with  Ms  = 11.1 
and  p0  = 5 torr,  radiation  cooling  is  not  negligible.  Nevertheless,  similar 
trends  of  the  flow  profiles  are  observed  in  their  end  wall  boundary  layers 
and  our  flat-plate  and  sidewall  boundary  layers. 


8-..  DISCUSSIONS  AND  CONCLUSIONS 

The  complete  set  of  partial-differential  equations  for  a laminar 
boundary  layer  in  ionizing  argon  have  been  solved  with  a six-point  implicit 
finite- difference  scheme.  The  new  features  in  the  analysis  are  the  inclusion 
of  the  radiation- energy  loss  and  the  appropriate  chemical  reactions.  The 
latter  also  include  the  atom- atom  reactions.  Account  was  taken  of  the  varia- 
tion across  the  boundary  layer  of  the  transport  properties  based  on  the  known 
elastic- scattering  cross-sections  for  an  argon  plasma.  The  compatibility 
conditions  and  the  electric- sheath  model  were  described  and  incorporated  into 
the  analysis.  The  flat  plate  and  shock- tube  sidewall  boundary- layer  flows 
were  analyzed  and  compared  with  interferometric  data  obtained  using  the  Iff  IAS 
10  cm  x 18  cm  Hypervelocity  Shock  Tube  equipped  with  a 23- cm  diam  Mach-Zehnder 
dual-wavelength  interferometer  at  shock  Mach  numbers  Mg  <*■* 13  and  16  at  an  ini- 
tial argon  pressure  pQ  ~5  torr  and  To  ~ 300  K. 

The  analysis  is  probably  the  most  complete  and  detailed  done  to 
date.  It  clearly  shows  that  the  measured  electron- numb er- densi ty  and  degree— 
of- ionization  profiles  in  the  boundary  layer  for  equilibrium,  frozen  and 
nonequilibrium  flows  are  more  sensitive  than  the  measured  complementary 
total-number  density  profile  for  determining  the  actual  state  of  the  boundary 
layer,  as  might  have  been  expected.  Agreement  between  theory  and  experiment 
for  the  density  profiles  appears  good  to  excellent  as  there  is  little  difference 
between  the  three  analytical  profiles  for  the  flat-plate  boundary- layer  flow. 
However,  agreement  of  experiment  with  analysis  for  electron-number  density  is 
only  fair  for  the  flat-plate  case.  The  experimental  data  lie  between  the 
analytical  frozen  and  nonequilibrium  profiles.  The  experimental  data  show  a 
bump  in  the  ne  profile  at  both  Mg  = 16.6  and  Mg  = 12.8  for  the  flat-plate 
cases,  while  analysis  only  predicts  a bump  for  the  frozen  case  at  Mg  = 12.8. 

The  nonstationary  character  of  the  inviscid  flow  may  have  an  important  effect 
on  the  boundary- layer  structure  at  high  shock  Mach  numbers  where  the  radiation 
energy  loss  is  significant.  However,  a substantive  explanation  for  this 
phenomenon  has  not  been  found. 

For  the  shock- tube  sidewall  boundary-layer  flow  it  was  shown  that 
there  is  no  bump  in  ne  profile  for  both  the  nonequilibrium  analysis  and  the 
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actual  experiments.  Better  agreement  with  experiment  was  found  in  this  case 
with  the  nonequilibrium  analysis,  rather  than  the  frozen  analysis.  A comparison 
of  the  indirect  temperature  measurements  and  theory  for  the  pj  asma-temperature 
profiles  shows  excellent  agreement. 

General  conclusions  obtained  from  the  comparisons  between  theoretical 
and  experimental  results  for  the  flat-plate  and  shock-tube  sidewall  boundary 
layer  flows  are  given  as  follows  (Liu  and  Glass,  Ref.  58): 

(a)  Near  the  wall  the  flows  are  in  nonequilibrium  or  near  frozen.  Equilibrium 
solutions  are  only  valid  for  the  flow  at  the  outer  part  of  the  boundary 
layer.  The  same  conclusion  was  drawn  by  KnDOs  for  a sidewall  boundary  layer 
and  by  Hutten  Mansfeld  for  a thermal  Rayleigh  boundary  layer. 

(b)  The  variation  of  the  flow  profiles  with  distance  x for  the  two  types  of 
boundary  layers  are  different.  For  example,  for  the  flat-plate  boundary 
layer  the  velocity  profile  is  almost  independent  on  x but  not  so  for  the 
sidewall  boundary  layer.  Consequently,  similarity  assumptions  are  reason- 
able for  the  velocity  and  atom  temperature  profiles  for  the  flat-plate 
case  where  the  freestream  variations  are  small;  and  even  for  the  sidewall 
velocity  profiles.  However,  large  errors  will  arise  in  electron-number- 
density  and  temperature  profiles  from  similarity  assumptions  in  both  cases. 

(c)  Unlike  the  total-number  density,  the  electron-nunber-density  profiles  are 
very  sensitive  indicators  for  comparing  frozen,  equilibrium  and  nonequili- 
brium analyses  with  experiments. 

(d)  The  thickness  of  the  electron  thermal  layer  is  of  the  same  order  of  magni- 
tude as  the  viscous  boundary  layer  for  both  flat-plate  and  sidewall  boundary 
layers . A similar  conclusion  was  made  by  Hutten  Mansfeld  from  a comparison 
between  his  analysis  and  Kuiper's  experimental  data  for  a thermal  Rayleigh 
boundary  layer. 
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FIG.  3 


Flow  at  t = tp  w 

SCHEMATIC  DIAGRAM  OF  A FLAT-PLATE  BOUNDARY-LAYER  FLOW  GENERATED  BY  A 
SHOCK  WAVE  (MOVING  FROM  LEFT  TO  RIGHT)  IN  A SHOCK  TUBE.  S - TRANSLATIONAL 
SHOCK  WAVE;  xe  - IONIZATION  RELAXATION  ZONE;  E - ELECTRON  CASCADE  FRONT; 

W - SHOCK-TUBE  WALLS;  Ex  - EXPANSION  WAVE;  B - BOUNDARY  LAYER;  M - FLAT- 
PLATE  MODEL;  R - RADIATIVE  COOLING  ZONE;  OTHER  QUANTITIES  AS  DEFINED  IN 
THE  TEXT. 
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VARIATION  OF  FREESTREAM  CONDITIONS  AT  EDGE  OF  BOUNDARY-LAYER  FLOW  WITH 
FLAT-PLATE  LEADING-EDGE  DISTANCE  x FOR  Mg  = 12.8,  pQ  = 5.01  TORR  AND 
T0  = 297  K. 
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5 INITIAL  FLOW  PROFILES  AT  x = 0 FOR  NORMALIZED  VELOCITY  f ' , DEGREE  OF 
IONIZATION  z,  ATOM  TEMPERATURE  0 AND  ELECTRON  TEMPERATURE  0 WITH  T) 
FOR  Mg  = 16.6,  p0  = 4.8l  TORR  AND  T0  = 296  K (CASE  1) . 


X = 0 


VARIATIONS  OF  TRANSPORT  PROPERTIES:  HEAVY  PARTICLE  PRANDTL  NUMBER  P 
RATIO  OF  DENSITY- VISCOSITY  PRODUCT  C,  SCHMIDT  NUMBER  Sc  AND  ELECTRON 
PRANDTL  NUMBER  Pre  WITH  T)  AT  x = 0 FOR  Ms  = l6.6,  pQ  = 4.8l  TORR  AND 
T0  = 296  K (CASE  1) . 
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FIG.  7 VARIATIONS  OF  NORMALIZED  FLOW  PROFILES  OF  VELOCITY  f',  DEGREE  OF 
IONIZATION  z,  ATOM  TEMPERATURE  9 AND  ELECTRON  TEMPERATURE  6 WITH 
LEADING-EDGE  DISTANCE  x FOR  Ms  = l6.6,  pQ  = 4.8l  TORE  AND  T0  = 2 
CASE  1) . 
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FIG.  8 VARIATION  OF  ELECTRON  PRANDTL  NUMBER  Pre  WITH  T)  FOR  LEADING-EDGE 
DISTANCES  x = 0 AND  x = 14  CM  FOR  Ms  = l6.6,  pD  = 4.8l  TORR  AND 
T0  = 296  K (CASE  1) . 
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FIG.  13  COMPARISON  OF  ANALYTICAL  AND  EXPERIMENTAL  PLASMA-DENSITY  PROFILE  p 
WITH  DISTANCE  y IN  THE  BOUNDARY  LAYER  FOR  Mg  = l6.6,  p0  = 4 .8l  TORR 
AND  To  = 296  K (CASE  l) . 


FIG.  14  COMPARISON  OF  ANALYTICAL  AND  EXPERIMENTAL  PROFILES  OF  NORMALIZED 

ELECTRON-NUMBER  DENSITY  rig/r^g  WITH  DISTANCE  y IN  THE  BOUNDARY  LAYER 
FOR  Ms  = 16.6,  p0  = 4.8l  TORR  AND  T0  = 296  K (CASE  l) . 
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COMPARISON  OF  ANALYTICAL  AND  EXPERIMENTAL  PROFILES  CF  NORMALIZED  DEGREE 
OF  IONIZATION  a/ctg  WITH  DISTANCE  y IN  THE  BOUNDARY  LAYER  FOR  Mg  = l6.6, 
Po  = **.8l  TORR  AND  T0  = 296  K (CASE  l)  . 
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FIG.  16  NONEQUILIBRIUM-FLOW  PROFILES  OF  NORMALIZED  VELOCITY  f ’ , DEGREE  OF 
IONIZATION  z,  ATOM  TEMPERATURE  9 AND  ELECTRON  TEMPERATURE  0 AS  A 
FUNCTION  OF  r\  AT  LEADING-EDGE  DISTANCE  x = l4  CM  FOR  Mg  = 12 .8, 
p0  = 5.01  TORR  AND  T0  = 297  K (CASE  2) . 
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FROZEN  FLOW  PROFILES  OF  NORMALIZED  VELOCITY  f',  DEGREE  OF  IONIZATION 
Z,  ATOM  TEMPERATURE  0 AND  ELECTRON  TEMPERATURE  0 AS  A FUNCTION  OF  T) 
AT  LEADING-EDGE  DISTANCE  x = l4  CM  FOR  Mg  = 12.8,  pQ  = 5.01  TQRR  AND 
T0  - 297  K (CASE  2) . 
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FIG.  18  COMPARISON  OF  ANALYTICAL  AND  EXPERIMENTAL  PLASMA-DENSITY  PROFILE  p 
WITH  DISTANCE  y IN  THE  BOUNDARY  LAYER  FOR  M_  = 12.8,  pQ  = 5. 01  TORR 
AND  T0  = 297  K (CASE  2) . 
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FIG.  19  COMPARISON  OF  ANALYTICAL  AND  EXPERIMENTAL  PROFILES  OF  NORMALIZED 

ELECTRON  NUMBER  DENSITY  ng/n^  WITH  DISTANCE  y IN  THE  BOUNDARY  LAYER 
FOR  Ms  = 12.8,  p0  = 5-01  TORE  AND  T0  = 297  K (CASE  2). 


FIG.  20  COMPARISON  OF  ANALYTICAL  AND  EXPERIMENTAL  PROFILES  OF  NORMALIZED 

DEGREE  OF  IONIZATION  cc/a&  WITH  DISTANCE  y IN  THE  BOUNDARY  LAYER  FOR 
Ms  = 12.8,  Po  = 5.01  TORR  AND  T0  = 297  K (CASE  2). 
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FIG.  21  SCHEMATIC  DIAGRAM  OF  A SHOCK-TUBE  SIDEWALL  BOUNDARY-LAYER  FLOW 
GENERATED  BY  A SHOCK  WAVE,  FIXED  AT  x = 0.  THE  WALL  IS  MOVING 
FROM  LEFT  TO  RIGHT  WITH  THE  SHOCK- WAVE  VELOCITY. 
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FIG.  23  NONEQUILIBRIUM- FLOW  PROFILES  OF  NORMALIZED  VELOCITY  F (OR  f ') , DEGREE 
OF  IONIZATION  z,  ATOM  TEMPERATURE  0 AND  ELECTRON  TEMPERATURE  0 AS  A 
FUNCTION  OF  y (OR  tj)  FOR  A SIDEWALL  BOUNDARY  LAYER  WITH  x = 9-5  CM, 

Mg  = 13.1,  p0  = 5.16  TORR  AND  T0  = 300  K. 


FIG.  24  COMPARISON  OF  ANALYTICAL  AND  EXPERIMENTAL  PROFILES  OF  PLASMA  DENSITY  P 
WITH  DISTANCE  y IN  THE  SIDEWALL  BOUNDARY  LAYER  FOR  Mg  = 13-1,  Po  = 5-l6 
TORR  AND  T0  = 300  K. 
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FIG.  25  COMPARISON  OF  ANALYTICAL  AND  EXPERIMENTAL  FROFILES  OF  ELECTRON-NUMBER 

DENSITY  ne  WITH  DISTANCE  y IN  THE  SIDEWALL  BOUNDARY  LAYER  FOR  Mg  = 13.1, 
p0  = 5-16  TORR  AND  Tn  = 300  K. 
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FIG.  26  COMPARISON  OF  ANALYTICAL  AND  EXPERIMENTAL  PROFILES  OF  DEGREE  OF  IONIZATION 
a WITH  DISTANCE  y IN  THE  SIDEWALL  BOUNDARY  LAYER  FOR  Mg  = 13.1,  p0  = 5-l6 
TORR  AND  To  = 300  K. 
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FIG.  27  BOUNDARY-LAYER  DISPLACE!®!®  THICKNESS  5*  AS  A FUNCTION  CF  DISTANCE  x FOR  THE  SIDEWALL  BOUNDARY 
LAYER  FOR  Me  = 13.1,  Po  = 5-l6  TORE  AND  T0  = 300  K. 
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DESCRIPTION  OF  COMPUTER  PROGRAM 


The  program  BLEIG  for  solving  Boundary  Layer  Equations  for  Ionizing 
Gases  is  written  on  CDC-6600.  The  other  versions  of  BLEIG  on  IBM- 370  and 
EDP-10  have  been  written.  The  CGS  unit  system  is  used  through  all  programs. 
The  flow  chart  of  the  program  is  given  in  Fig.  28.  Free  format  of  the  input 
data  is  used.  The  computer  program  is  given  in  Appendix  F. 

A.l  Main  Program 

Before  describing  the  notations  used  in  the  main  program,  the  main 
features  of  the  calculation  procedure  will  be  reviewed.  A rectangular  grid 
system  indicated  in  Fig.  1 has  been  adopted,  the  j-lines  running  in  the  i) 

(or  y)  direction,  i.e.,  normal  to  the  plate,  and  the  i-lines  in  the  | (or  x) 
direction,  i.e.,  parallel  to  the  plate.  Conditions  along  seme  initial  j-line 
are  known  and  the  conditions  along  the  (j+l)-line  have  to  be  determined.  The 
main  steps  in  the  procedure  are  then, 

(1)  If  this  is  an  original  run,  the  solutions  of  the  flow  at  leading  edge 
(x=0)  are  obtained  by  calling  SUBROUTINES  SETUP  and  BEGIN.  If  this  is 
not  the  original  run,  then  the  solutions  at  some  particular  point  of  x 
are  read  in  from  a restart  file  in  a magnetic  tape. 

(2)  From  the  known  solutions  of  F,  z,  0 and  9 on  the  j-line,  the  solutions 
at  (3+1) -line  are  assumed. 

(3)  Using  assumed  solutions  at  ( j+l)-line,  the  new  velocity  profile  F at 
(j+l)-line  is  calculated. 

(4)  Using  new  value  F and  assumed  z,  0 and  9 values,  the  new  value  of  z is 
calculated. 

(5)  Using  new  values  of  F and  z and  assumed  values  of  0 and  9,  new  value  of 
0 is  calculated. 

(6)  Using  new  values  of  F,  z and  0 and  assumed  9 value,  new  9 value  is 
calculated. 

(7)  Repeat  from  step  (3)  to  step  (6)  until  the  solutions  at  (j+l)-line 
converge  to  satisfy  a preset  criterion. 

(8)  Calculate  the  boundary-layer  characters  and  determine  a suitable  step- 
size  Ax. 

(9)  Use  the  same  procedure  to  advance  from  the  ( j+l)-line  to  the  (j+2)-line 
and  so  on. 

The  notations  used  in  the  main  program  are  listed  and  explained 

below: 


A-l 


EDU  (MREAD) 

= 

Ug  values  at  the  edge  of  boundary  layer. 

EDTA  (MREAD) 

= 

Tag  values  at  the  edge  of  boundary  layer.. 

EDTE  (MREAD) 

= 

Teg  values  at  the  edge  of  boundary  layer. 

EDALP  (MREAD) 

= 

dig  values  at  the  edge  of  boundary  layer. 

EDP  (MREAD) 

= 

Pg  values  at  the  edge  of  boundary  layer. 

XDIS  (MREAD) 

= 

Reading  x values  for  the  freestream  conditions. 

DENNE  (I) 

= 

ng  at  grid  point  (i,  j+l) . 

WF  (I) 

= 

F at  grid  point  (i,  j+l) . 

WZ  (I) 

= 

z at  grid  point  (i,  j+l) . 

WTA  (I) 

= 

9 at  grid  point  (i,  j+l) . 

WTE  (I) 

= 

9 at  grid  point  (i,  j+l). 

SF  (I) 

= 

f at  grid  point  (i,  j+l) . 

SFX  (I) 

= 

f^  at  grid  point  (i,  j+l) . 

WFP  (I) 

= 

F at  grid  point  (i,  j) 

WZP  (I) 

= 

z at  grid  point  (i,  j) 

WTAP  (I) 

= 

9 at  grid  point  (i,  j) 

WTEP  (I) 

= 

9 at  grid  point  (i,  j) 

r*  r 

GUESS  (4) 

~ 

Initial  guess  values  of  [Cf']w,  z'  J , 

and  9 • 
w 

WORK  (I) 

= 

A function  evaluated  at  (i,  j+l). 

f(i) 

= 

Weighted  F value  at  (i,  j+l) . 

Z (I) 

= 

Weighted  z value  at  (i,  j+l) . 

THETA  (I) 

= 

Weighted  9 value  at  (i,  j+l). 

THETE  (I) 

= 

Weighted  9 value  at  (i,  j+l) . 

STAY  (I) 

= 

tj  at  (i,  j+l) . 

YREAL  (I) 

= 

y at  (i,  j+l) . 

CR  (I) 

= 

C at  (i,  j+l). 

A- 2 


CKPR  (I) 

= 

C/Pr  at  (i,  j+1) 

CRSC  (I) 

= 

C/Sc  at  (i,  j+,1) 

CRPRE  (I) 

= 

C/Pre  at  (i,  j+1) 

RHO  (I) 

= 

p at  (i,  j+1) 

CRP  (I) 

= 

C at  (i,  j) 

GRPRP  (I) 

= . 

C/Pr  at  (i,  j) 

CRSCP  (I) 

= 

C/Sc  at  (i,  j) 

CRPREP  (I) 

= 

C/Pre  at  (i,  j) 

OLDWF  (I) 

= 

Previous  iterative  value  of  F at  (i, 

J+1) 

OLDWZ  (I) 

= 

Previous  iterative  value  of  z at  (i, 

j+1) 

OLDWTA  (I) 

= 

Previous  iterative  value  of  0 at  (i, 

j+D 

OLDWIE  (I) 

= 

Previous  iterative  value  of  0 at  (i, 

j+1) 

LAMD 

= 

* = Ua  + \) 

LAMDE 

= 

\ 

MJ 

= 

n 

K 

= 

k 

HA 

= 

na 

HE 

= 

ne 

MOW 

= 

HEDOTA 

= 

(ne)a/ne 

HEDOTE 

= 

(ne)e/ne 

UED 

= 

U5 

PED 

= 

P5 

TAED 

= 

Ta6 

TEED 

= 

Te6 

ALPED 

= 

a5 

TW 

= 

Tw 

A-3 


zw 

= 

z 

w 

uw 

= 

uw 

ROLMLJE 

= 

P6P6 

CN 

= 

A 

DETA1 

• = 

N 

= 

Maximum  value  of  i 

DEX 

= 

Mi 

IBC 

— 

If  IBC  = 1,  sheath-wall 

If  IBC  - 2,  sheath-wall 

A(I) 

= 

Matrix  element  A^ 

B(I) 

= 

Matrix  element  B^ 

C(I) 

= 

Matrix  element  Ci 

D(I) 

= 

Matrix  element  Di 

D1(I) 

= 

l/fd+k^k1] 

D2(I) 

= 

2/[(l+k)k2i"1A*i12] 

X(t) 

= 

Xi 

MA 

= 

m 

a 

R 

= 

R 

TION 

= 

TI 

TEXC 

= 

T* 

SAA 

= 

S* 

aa 

SAE 

= 

S* 

ae 

TACRI 

= 

T 

a,cn 

TBCRI 

= 

T 

e jCri 

QAAC 

= 

Used  in  o-  = QAAC/T°'25 
8.8.  8. 

QAIC 

= 

Used  in  a = QAIC/T0,09 
8.1  & 

is  used 
and  a are  used 


A-1+ 


CC1,  CC2,  CC3, 
CC4,  TEEM.,  DD1, 
DD2,  DD3,  DI>4  = 


EEL,  EE2,  EE3 


HVC 

ZEET 

DX 

XMAX 

EPS 

ITYPE 

ICASE 

CONTL 

YMM 

THICKD 

J 

XREAL 

XI 

ITER 

IFRNTX 


Used  in  the  equation: 

o-  = (CC1  + CC2-T  + CC3‘T  2 + CC4-T  3)  x 10_l6 

08,  G 6 G 

if  T < TEEA 

G — 

ct  = (EDI  + DD2-T  + DD3-T  2 + DD4-T  3)  x 10~l6 

08,  G G G 

if  T > TEEA 
e 

Used  to  calculate  the  statistical  weight  of  ion,  Zj, 
see  FUNCTION  EQK:  ZI  = EE1  + E£2/exp(EE3/T) 

vc 

Zeff 

Ax 

Maximum  value  of  x 
Tolerance  criterion 

ITYPE  = 1 for  flat-plate  boundary  layer 
ITYPE  = 2 for  sidewall  boundary  layer 
ICASE  = 1 for  nonequilibrium  flow 
ICASE  = 2 for  frozen  flow 

Parameter  used  in  the  control  of  step  size  Ac 
y value  in  mm 

displacement  thickness  5*  in  cm 
J 

x in  cm 

I 

Iteration  number 

Solution  step  print  selector  in  x direction.  Output 
will  be  printed  on  the  line  printer  every  IHOTX  step 
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IPRNTY 


UTLE 


Solution  step  print  selector  in  y direction.  Output  will 
be  printed  on  the  line  printer  every  IPRNTY  step 
Restart  file  solution  step  selector.  A restart  record 
will  be  written  to  the  restart  file  every  IFILE  step 


FIN 

= 

Read  in  restart  file  name 

FONT 

= 

Written  out  restart  file  name 

MREAD 

= 

Number  of  points  freestream  conditions 

N1 

= 

N-l 

DWFDY 

= 

(3F/dt))w 

DWZDY 

= 

(dZ/dti)w 

DWTADY 

= 

(de/dr))w 

DWTEDY 

= 

( d0/dr|  )w 

GAML 

= 

p^6U5  at  (N,  J) 

GAM2 

= 

P5^6U5  at  ^N’  d+1^ 

TAUX 

= 

Te6//Ta6  = T 

AL1AU 

= 

T 

0 

TAUU 

= 

BETAF 

= 

Pf 

BETAZ 

= 

BE1ATA 

= 

^a 

BET  ATE 

= 

^e 

TEW 

= 

^ew 

sew 

= 

Scw 

SLOPE 

= 

z'  or  0' 
w w 

l'A 

= 

T 

a 

TE 

= 

Te 

XIUED 

= 

*£1/  Ug 
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A.  2 Subroutines 


Subroutine  IOTERP: 

The  Aitken-Lagrange  and  Lagrange  interpolations  have  been  applied 
before.  However,  it  was  found  that  these  interpolations  are  not  suitable  for 
the  sidewall  boundary  layer  where  the  variation  of  a§  with  x is  significant. 
Therefore,  a linear  interpolation  IOTERP  is  applied  in  the  present  program. 

Y = The  resulting  interpolated  function  value. 

X = The  argument  value  specified  by  input. 

M = An  input  value  which  specifies  the  number  of  points  in  table  (XF,  YF) . 

XF  = The  input  vector  of  argument. 

YF  = The  input  vector  of  function  values  of  the  table. 

Subroutine  SETUP: 

Equations  (60)  - (63)  are  the  ordinary  differential  equations  with 
two-point  boundary  conditions.  The  Newton-Raphson  method  for  the  iteration 
techniques  (see  Subroutine  BEGIN)  is  applied.  In  the  BEGIN  subroutine,  four 
guess  values  for  [Cf"]w,  [(C/Sc)z']w,  [(C/Pr)e']w  and  6w  (or  [Cf"]w,  Z*, 
[(C/Pr)e']w  and.  0^)  are  required.  If  the  guess  initial  values  are  not  reasonable, 
then  the  commutation  time  is  more  consumable  and  stable  solutions  cannot  be 
found  in  the  Subroutine  BEGIN.  In  order  to  avoid  this  difficulty,  same  knowledge 
about  how  to  guess  the  initial  values  are  necessary. 

From  the  integral  method  described  in  Ref.  10,  the  following  ways  to 
guess  the  initial  values  are  presented: 


GUESS  (1)  = Cw  fj  (A1) 

where 

f;  - vat5 

- 2(1  - 7) 

y - Vu6 


\ - 2“l/(Ff5) 


p “’ll  ac  a.  a*  PIT 

G = ai  + — + -9-  + IT-  + — + 3 V5  + 3 *1*6  + Tf  a4*5 


+ I a4a6  + 5 a5a6 


a-L  = 2(1  - 7) 

% = -5(1  - 7) 


= 6(1  - 7) 


a6  = -2(1  - 7) 


guess  (4)  = ft 


Using  the  boundary  condition  0^  = 0 and  Eq.  (63),  we  obtain  GUESS  (4)  = 1.0. 

In  the  above  equations,  the  transport  properties,  C,  Sc,  Pr  are 
evaluated  at  some  particular  distance  t)*  from  the  wall.  However,  it  is 
difficult  to' determine  tj*  from  theory.  From  our  experience,  the  guess 
transport  properties,  are  approximately  equal  to  the  average  values  of 
that  at  the  f reestream  and  at  the  wall . 


Subroutine  INERG : 

This  is  a subroutine  to  perform  the  following  integration  by  three- 
point  difference  formulae . ' 


-J"  — 


i=0 


L + k)  y ^ + ^ 


It  is  worth  noting  that  I = 1 in  the  program  corresponds  to  i = 0 
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Subroutines  BEGIN  and  BLF: 


These  subroutines  are  used  to  solve  Eqs.  (60)  - (63)  by  means  of 
the  Newton-Raphson  iteration  techniques  (Ref.  10) . The  following  notations 
are  applied: 


NORD  = Number  of  ordinary  differential  equations 


i » 1 to  9 


Y (8)  = e 

Y (9)  - ~ 9’ 

e 

YFRIM  (i)  = Y'(i) 

Subroutine  RKGIL: 

This  subroutine  uses  the  Runge-Kutta  method  with  Gill's  coefficients 
for  the  solution  of  initial- value  problems. 

Y = Input  vector  of  initial  values,  y 

YPRIM  = y' 

NORD  = Nunber  of  equations 

A = Lower  bound  of  the  interval 

B = Upper  bound  of  the  interval 

STEP  = Step  size 

DER  = An  auxiliary  storage  array  of  y ' 


Subroutine  DENS  IT: 

This  subroutine  is  used  to  calculate  na, 

Te  and  a. 


n6  and  p by  given  p,  Ta, 


Subroutine  TRANSP: 

This  subroutine  is  used  to  calculate  the  transport  parameters . The 
following  notations  are  used: 


QAA 

= 

cr  „ 

aa 

QAI 

zz 

cr  . 

ai 

QEA 

— 

ea 

QEI 

= 

aei 

QEE 

— 

cr 

ee 

QII 

= 

aii 
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CA 

CE 

MU 

LAMD 

LAMDE 

DAW 

VEAI 


_ r ^5Tt 


i 


m 


f 


h 

\ 

\ 

D 


amb 


.a  + n or  . ) 
a ea  e ei' 


Subroutine  RAZES: 

This  is  a subroutine  to  calculate  the  reaction  rates. 


NEDOTA 

= 

(ne)a/ne 

NEDOTE 

= 

(n  ) /n 
v e'  e'  e 

KFA 

= 

kfa 

KRA 

= 

kra 

EQK 

= 

Keq. 

KFE 

= 

kfe 

KRE 

= 

kre 

SAA 

— 

First  excitational 
atom  and  atom 

SAE 

— 

First  excitational 

atom  and  electron 
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Function  EQK: 


This  is  a function  to  calculate  the  equilibrium  rate  constant  for 
a given  temperature . 

ZI  = Statistical  weight  for  an  ion 


Function  QRAD: 

This  function  is  used  to  calculate  the  radiation  energy  loss. 

HVC  = Cut-off  frequency  of  the  plasma 

ZEFF  = Effective  nuclear  charge 

QC  = Continuum  radiation  energy  loss 

QEAD  = Total  radiation  energy  loss 

NE  = Electron  number  density 

TE  = Electron  temperature 


Subroutines  MINVRS  and  SUBJCS: 

These  subroutines  are  used  to  calculate  the  matrix  inverse  of  matrix 
A.  The  input  matrix  and  restating  matrix  inverse  are  specified  by  A. 

Subroutines  CRANK  and  TRDG: 

Subroutines  CRANK  and  TRDG  are  an  algorithm  to  solve  a tridiagonal 
matrix  most  efficiently.  The  matrix  equation  is  given  by 

AiWi-l  + BiWi  + CiWi+l  = Di 

where  i = 1 to  N-l,  and  WQ  and  Wjj  are  given.  The  matrix  form  is 


BI  C1 
*2  B2  °2 
A3  B3  °3 


S-2  BN-2  CN-2  WN-2 

^-1  BN-1  WN-1 


D£  = Di  - A^ 


If  (Sw/&n)Q  is  given,  then 

B!  - \ + Bi 


For  the  present  boundary  layer  case,  WK  (=1)  is  given,  and 
therefore,  - C^. 

The  arguments  of  CRANK  are 

IBC 

= Type  of  boundary  condition  at  the  wall;  1 for  Dirichlet 
type  boundary  condition  and  2 for  Neumann  type  boundary 
conditions 

SLOPE 

n 

si* 

o'"  y 

DETA1 

II 

* 

WSTART 

= Resulting  solutions 

Subroutine  COEFF: 

This  is  the  subroutine  used  to  calculate  the  matrix  element  A*, 
Bi,  Ci  and  Di  (given  by  Eqs.  75)  used  in  Subroutine  CRANK. 

Subroutine  DXCTL: 

This  is  the  subroutine  used  to  control  the  step  size  ^c.  The 
argument s are 

DX  = Step  size  Zfic 


THICKD  = Displacement  thickness  of  boundary  layer  5* 


DIET  = Maximum  value  of  |w(i,  j+l)  - W(i,  j)  | 

HER  = Iteration  number  used  for  previous  step  size 

The  step  size  DX  is  given  by 

DX  = CONTL  • |5*  | 

Subroutine  OUTDSK:  Subroutine  to  write  out  the  results  on  Tape  3 • 

The  subroutine  OUTDSK  is  used  for  writing  the  necessary  results  on 
the  magnetic  tape  under  the  output  restart  file  name.  These  results  are 
needed  for  the  calculation  starting  at  x greater  than  zero. 

A. 3 Input  and  Output 

The  input  data  for  BLEIG  is  a description  of  boundary  layer  charac- 
ters, freestream  conditions,  physical  parameters  and  parameters  of  boundary 
layer  structure.  Where  possible,  'free  format'  input  has  been  used  to  sinplify 
the  use  of  default  values. 

The  following  cards  are  needed  for  input  data  (using  CGS  units): 

(1)  ITYPE,  ICASE,  IBC,  IPRNTX,  IFRNTY,  IF HE 

(2)  N,  DETA1,  XREAL,  DX,  XMAX,  K,  CN,  TW,  EPS,  CONTL,  ZW 

(3)  MREAD — NUMBER  OF  POINT  FREESTREAM  CONDITIONS  READ 

(4)  XDIS (MREAD),  EDALP ( MREAD ) , EDU(MREAD),  EDTE(MREAD),  EDTA(MREAD), 

EDP(MREAD) 

(5)  fin,  FOUT  — READ  IN  AND  READ  OUT  RESTART  FILE  NAMES 

(6)  ISAS  —IF  IGAS  EQUALS  TO  ZERO,  THEN  DEFAULTED  VALUES  ARE  USED.  IF  IGAS 
IS  NOT  EQUAL  TO  ZERO,  THEN  READ  THE  FOLLOWING  DATA 

(7)  MA,  R,  TION,  TEXC,  SAA,  SAE,  TACRI,  TECRI 

(8)  QAAC,  QAIC,  CC1,  CC2,  CC3,  CC4,  DD1,  DD2,  DD3,  DD4,  TEEA 

(9)  EE1,  EE2,  EE3 

(10)  HVC,  ZEFF 

For  example,  the  following  six  data  cards  are  used  for  a nonequilibrium 
flat-plate  boundary  with  Ms  = 16.6,  p0  = 4.8l  torr  and  T0  = 296  K: 

(1)  1,  1,  1,  5,  1,  200 

(2)  70,  0.035,  o,  l.E-6,  14,  1.05,  0.75,  296,  1.0E-4,  0.1,  l.E-5 

(3)  1 

(4)  14,  0.021,  4.86E5,  1.049E4,  1.049E4,  0.27E7 


The  following  outputs  are  given  in  the  present  program: 


(6)  If  XREAL  in  input  data  is  zero,  then  the  initial  GUESS  value  obtained 
from  SETUP  and  Y(l)  to  Y(9)  obtained  from  BEGIN  are  given. 

(7)  XREAL,  XI,  J,  ITER,  DX,  UED,  TAED,  TEEP,  ALP ED  and  PED 

(8)  I,  ETAY,  Y(nm),  WF,  W2,  VJTA,  WTE,  CR,  CRPR,  CRSC,  CRPKE,  RHO  and  NE 


I 
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IS  ITER  > 15? 

YES  1 YES 

CALCULATE  BOUNDARY- LAYER  CHARACTERS  « wAPMTNr  ' IS  IT  ACCEPTABLE? 


SELECT  OUTPUT, RESTART  FILE  AND  STEP-SIZE Dx 


IS  x > x ? 
max 


FIG.  A.l  FLOW  CHART  OF  THE  PROGRAM  BLEIG. 


APPENDIX  B 


THERMODYNAMIC  QUANTITIES 


The  definitions  and  the  equations  for  thermodynamic  quantities 
of  a two-temperature  ionizing  monatomic  gas  will  be  summarized.  The  detailed 
formulations  of  the  thermodynamic  equations  can  be  found  in  the  text  books  and 
in  the  literature. 

In  this  study,  we  are  dealing  with  a mixture  of  atoms  (a),  ions  (i) 
and  electrons  (e) . Atoms  and  ions  have  the  same  teniperature  Ta  and  electrons 
have  a temperature  Te  For  each  species  j we  have 


’ "d  TJ 

Tj 


V/ 


°P3  + h3 


where  p-j  is  the  pressure  of  j species,  n^  is  the  nunfcer  density  of  j species, 
kg  is  tne  Boltzmann  constant,  hj  is  the  enthalpy  per  unit  mass  of  j species, 
Cpj  is  the  specific  heat  at  constant  pressure  per  unit  mass  of  species  j and 

hj  is  the  chemical  enthalpy  per  unit  mass  of  species  j . 

The  total  pressure  of  the  mixture  is  defined  by  Dalton's  law  as 


and  the  density  of  the  mixture  p is  given  by 


where  p 
tron  is' 


pj  = m-jnj,  DLj  is  the  mass  of  the  species  j.  Since  the  mass  of  an  elec- 
s negligibly0  small  compared  with  the  mass  of  an  atom,  p is  written  as 


p - ma(»a  * ne) 


where  = ng  for  singly  ionized  gas . 

The  mass  fraction  of  ions,  a,  which  is  also  called  the  degree  of 
ionization,  is  given  by 


p n + n 
K a e 


B-l 


r 


I 


The  mass  fraction  of  electrons  is  given  by 


o in 
— — e ct 


P m 


and  the  mass  fraction  of  atcms  is 


— = 1 - a 

P 


The  total  number  density  of  the  plasma  is 
n = ) n . = n + 2n 

L J a e 

n 

The  relations  between  n , n and  a are  given  by 
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n = p(l  -a)/] 


m 


ne  = P°^ma 


Therefore,  n can  be  written  as 


n = p(l  + a)/m 


The  total  pressure  p of  the  plasma  becomes 


5 " (na  + ne)kB  Ta  + ne  S T< 


= pR(Ta  + QTe) 


(B7) 


(b8) 


(B9) 


(BIO) 

(Bll) 


(B12) 


(B13) 


where  R is  the  gas  constant  referred  to  the  atomic  gas  and  defined  by  R = 

Let  Nj  be  the  number  of  j -particles  in  the  volume  V and  F.  be  the 
Helmholtz  free  energy  function  of  one  particle  of  the  species  j,  then  the 
Helmholtz  free  energy  F of  the  mixture  is  given  by 


* 


F 

J d 


(BlU) 
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where  Nj  = n^V  and  F^  is  given  by 

Fi  = "V^168  gJ  ' lo8  Nj  + lo6  + 1] 


(B15) 


t . i 

gd  80(1  gj  8X6  the  partition  functions  associated  with  the  translational  and 
internal  degree  of  freedom,  respectively. 


The  entropy  S is  defined  by 


A 


where  Sj  = -(SF^/StD^ )y. 

The  internal  energy  E is  defined  by 


E = F 


+ TS  = Y H 


dEd 


where 


where 


The  enthalpy  H is  defined  by 

9 -XBd  sj 

A 


Hd  = Ed  + pd  % 


Finally,  the  Gibb's  free  energy  of  the  plasma  G is  given  by 

0=XVj 


A 


where 


(Bl6) 


(bit) 


(B18) 


(B19) 


°4  - 94  - TSJ 


According  to  the  above  relationships,  the  specific  properties,  i.e., 
the  properties  per  unit  mass  of  an  i deal  ionizing  monatomic  gas  can  be  derived 
from  the  partition  functions  of  each  species  involved.  The  partition  function 


associated  with  the  translational  mode  is  given  by 


gd  = 


(2um.iT.) 

— 


3/2 


6oj 


'^3^3 


(B20) 


where  €0.  is  the  excess  energy  of  the  ground  states  of  the  species  above  the 
reference  (ground)  energy  level,  gD.  is  the  probability  or  the  statistical 
weight  of  the  ground  energy  level  J eoj  and  h is  the  Planck  constant. 

The  internal  mode  of  electronic  excitation  is  always  assumed  to  be 
at  its  ground  state.  The  internal  partition  function  for  electronic  excitation 
is  given  by 


gJ 


(B21) 


n 


where  en  is  the  energy  of  the  n-th  state  of  the  particular  species  above  its 
ground  state  and  gn  is  the  statistical  weight  of  state  n. 

Using  Eqs.  (B20)  and  (B21) , the  specific  internal  energy  e (e  = E/pV) 
and  the  specific  enthalpy  h (h  = H/pV)  of  the  mixture  can  be  derived  as  follows: 


= | R(T  + QTj  + GfflT_ 


h - e + £ = | R(T  + OT  ) + QRTt 

n c.  a G ± 


(B22) 

(B23) 


The  component  specific  heats  CL.  and  the  frozen  specific  heat  of  the 
mixture  Cpf  at  constant  pressure  are  defined  by 


‘I 


"pd 


"pf 


J 'V 


i c 


(B24) 


pd 


given  by 


The  specific  heats  are  all  5%/2  per  particle.  Therefore,  Cpf  is 


V ' I E(1  * 


(B25) 


Similarly,  the  specific  heats  and  the  frozen  specific  heat  of 
the  mixture  Cyf  at  constant  volume  are  defined  by 


■(fti 


(B26) 

Contd. 
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uvf 


R(1  + a) 


(B26) 


When  the  mixture  is  in  chemical  equilibrium,  Gibb ' s free  energy  has 
its  minimum  for  all  possible  changes  in  composition  of  a system  at  a given 
pressure  and  a given  temperature . The  equilibrium  equation  for  the  degree  of 
ionization  auq  can  be  obtained  by  using  this  condition.  A detailed  derivation 
of  aeq  = °eq(P»  T)  Based  on  this  concept  was  given  by  Glass  and  Takano.  Another 
derivation  of  the  Saha  equation,  or  equilibrium  equation,  is  described  as 
follows:  Let  the  overall  reaction  paths  be  represented  by 


A + A 


A + e 


-ft 


ra 


fe 


re 


A + e + A 


A + e + e 


(B27) 


where  A denotes  atoms. 

Production  rates  of  electron  number  density  associated  with  Eq. 
(B27)  can  be  expressed  as  follows: 


(rijn  = k-  n 2 - k n n c 
e'a  fa  a ra  a e 


(B28) 


(n  ) = k n n - k n ~ 

' e'e  fe  a e re  e 


(B29) 


where  kf  and  kr  are  forward  and  reverse  rate  coefficients,  respectively. 

In  an  equilibrium  state,  the  forward  and  reverse  reactions  are  balanced 
in  the  process.  Therefore, 


*fa(T  ) 

- VV  ■ i^rry 

ra  a'  ^ a,eqv  a' 


n 


(B30) 


kfe(T)  ne  ea(Te) 

*3*7  ■ V V ■ 

where  K^T)  is  the  equilibrium  constant,  which  is  given  by 


(B31) 


2g  / 2 m k T n3/2 

VT)  ■ ef  ( — r6-  ) «P  (-TXA) 


(B32) 
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and  6a>  8t  are  the  electronic  excitation  partition  functions  of  atom  A 
and  ion  A*  particles. 

Using  Eqs.  (BIO)  and  (Bll) , the  equilibrium  value  of  degree  of 
ionization,  aeq,  is  obtained  from  Eqs.  (L30)  and  (B31)  as  follows: 


a m 

rrV-  - 7T  WT) 

eq  peq  H 


(B33) 


where  T stands  for  T and  T . 

a e 

If  the  system  is  in  thermal  equilibrium  (Te  = Ta  = T),  then  the 
equilibrium  value  of  density  pgq  is  given  by 


peq  ~ RT^f 


Ta  7 

eq' 


(B34) 


Introducing  a characteristic  density  for  ionization,  pj,  the 
equilibrium  equation  for  a becomes 


aeq(p>  T) 


p / TI  \5/2  Tj/T 

) e 


(B35) 


where 


pi 


= 2 


“e  ^/2  Sj  (27^Ti)3/2 
my  g 73 

a / 6r  h 


m5/2 


(B36) 


For  argon,  pj  = 150.27  g/cm3. 

Equation  (B35)  is  known  as  Saha's  equation.  Given  pressure  p and 
temperature  T,  the  degree  of  ionization  a is  calculated  from  Eq.  (B35)  if  the 
flow  is  in  chemical  equilibrium. 


Production  rates  of  electron  number  density,  (ne)a  and  (h»)e,  given 
by  Eqs.  (B28)  and  (B29)  became 


<4e>a  ' kra<Ta>  ".tVT.)a.  ' <B37) 

<Ve  ‘ kre<Te>  WTe>Ha  ‘ n»1  G>38) 

The  equilibrium  degree  of  ionization  Q£e  is  obtained  by  setting  the 
square  brackets  in  Eqs.  (B37)  and  (B38)  equal  to  zero  and  solving  for  a.  Thus 
the  equilibrium  degree  of  ionization  is  a function  of  the  local  temperature 
and  pressure.  This  does  not  imply  that  (ne)a  = 0 and  (rie)e  = 0 or  even  that 
ne  approaches  zero  at  equilibrium,  since  to  achieve  equilibrium  the  reverse 
reaction  rate  coefficient  kra  and  kre  must  become  infinite.  As  kra  and  kre 
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approach  infinity,  equilibrium  is  attained.  The  case  k_H  = 0 and  k , . « o 
corresponds  to  frozen  flow.  The  actual  value  of  (6Je7as  a function  of  the 
boundary  layer  coordinate  tj  cannot  be  calculated  frcmfoe  chemical  kinetics 
but  must  be  determined  by  solving  the  boundary  layer  equations. 

In  the  inviscid  flow  region  if  the  flow  is  frozen  then  a is  constant 
since  the  antoipolar  diffusion  velocity  of  ions  and  electrons  is  assumed  very 
small  and  can  be  neglected.  However,  if  the  flow  is  frozen  then  a is  not 
constant  in  the  boundary  layer  flow  due  to  the  diffusion  process  of  ions  and 
electrons. 

The  following  mathematical  expressions  for  the  flow  conditions  are 

given: 


Equilibrium: 


(nj 


lim 


e'eq  k -*«c  [(Ae^a  + ^ne^ 

Kt 


a -»a 


eq 


Frozen: 


frozen  " 0 


Nonequilibrium: 


APPENDIX  C 


With 


i = 1,  W^  = F: 


With  i = 2,  W^  = z: 


With  1=3,  W^  ° e: 


EXPRESSIONS  OF  FUNCTIONS  x 


U) 


i - 


= c 


*2  - C + f + 2|  f 


x3  = -0f  F 


x,  = 2J  F 


-r  (?) 


1 Sc 


■ft). 


+ f + 2|  f, 


x_  = “0  F + T 
3 z a ne 


x2j  = 2|  F 


x5  = 0 


X1  Pr 


+ f + 2fc  f . 
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APPENDIX  D 


EFFECTS  OF  CHEMICAL  REACTIONS  ON  BOUNDARY  LAYER 


The  large  change  of  the  chemical  reaction  rates  with  temperatures 
has  an  important  effect  on  the  boundary  layer  structure  for  the  case  of  large 
degree  of  ionization.  Two  important  effects,  mathematical  and  physical  effects, 
of  chemical  reactions  on  the  boundary  layer  structure  are  described  as  follows: 

1.  The  equation  for  conservation  of  electrons  is  given  by  Eq.  (21) . 

This  nonlinear  equation  is  coupled  to  the  equations  of  conservation  of  plasma 
and  electron  energies,  momentum  and  mass  of  the  remaining  species  through  the 
functions  C,  Sc,  f,  f'  and  . Without  any  loss  in  the  general  character  of 
the  equation  for  electron  species, we  assume  a local  similarity  equation  with 
C/Sc  = 1,  f = 1 and  f ' =0: 


where 


z"  + z*  + K'  n =0 
e 


m 


(Dl) 


K'  = -21 

2 

Ws 


A?  p“6 


In  order  to  illustrate  the  nature  of  the  problem  involved,  the  following 
linearized  and  simplified  version  of  Eq.  (Dl)  is  considered: 


z"  + z'  - K z = 0 


where  K is  a const  suit  which  includes  the  chemical  reaction  term. 


(D£) 


This  simplified  and  linearized  equation,  Eq.  (D2) , contains  till  the 
difficulties  associated  with  the  integration  scheme  used  in  numericsil  solution 
of  the  boundary  layer  equations.  The  general  solution  of  Eq.  (D2)  is 

1 + \/l  + 4K  . n/1  + 4k  - 1 
- £ T|  g ’l 


z = C1  e 


+ C2  e 


(S3) 


where  K>  0.  The  constants  Ci  and  C2  can  be  determined  from  the  two-point 
boundary  values  of  z = zw  at  t)  = 0 and  z - 1 at  i)  -too  or  z1  = zy  at  q = 0 and 
z = 1 at  t|  «♦  oo. 

It  is  seen  from  Eq.  (d3)  that  if  the  initial  guess  for  is  not 
correct,  then  the  second  term  in  Eq.  (D3)  will  dominate  at  large  values  of  t| 
so  that  z -» oo  at  t)  ->  oo.  Only  if  K = 0,  z at  ij  -+»  will  achieve  a finite  vsilue 
that  can  be  used  to  refine  the  guess  for  z^. 


D-l 
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Blottner  (Ref.  17)  has  used  a finite- difference  method  to  solve 
Eq.  (21) . However,  he  did  not  linearize  the  source  term  he  as  we  have  done  in 
Chapter  5.  The  solutions  given  by  Blottner  (Ref.  17)  appear  to  be  cases  for 
which  K « 1.  The  method  used  for  linearization  in  Chapter  5 is  similar  to  the 
method  of  quasilinearization  used  by  Fay  and  Kaye  (Ref.  5)  in  the  solution  of 
similar  nonequilibrium  boundary  layers . 

For  the  frozen  flow  (he  = 0) , Eq.  Cl  becomes: 


z"  + z'  = 0 


(04) 


Solution  of  Eq.  04  is  given  by 


z = C1e"1fJ  + c2 

Therefore,  no  mathematical  difficulty  associates  with  the  frozen  flow. 


2*  The  dependence  of  the  reaction  rate  coefficients  on  the  temperatures 

are  described  in  Section  3.3.  For  low  temperature  (below  1000  K) , kfa(Ta) 
and  kfe(Te)  are  very  small  collared  with  kra(Ta)  and  kre(Te),  respectively, 
and  the  forward  reaction  rates  can  be  neglected.  With  temperatures  above 
10,000  K,  forward  reaction  rates  become  significant.  However,  these  forward 
reaction  rates  are  still  very  small  compared  with  the  reverse  recombination 
rates.  The  following  relations  are  satisfied  at  temperature  about  10,000  K: 


kre  V » kfe  \ % 


2 2 
k n n » k n 
ra  a e fa  a 


■5  p 

k n J > k nn 
re  e ^ ra  a e 


Therefore  the  reverse  recombination  rate  due  to  electron-ion-electron  collisions 
is  the  important  process  in  the  boundary  layer  structure. 

However,  in  the  region  where  atom  temperature  is  about  25,000  K (for 
example,  in  the  region  near  the  shock  front  where  atom  temperature  is  much 
greater  than  electron  temperature) , the  following  relations  hold: 


2 2 
k_  n » k n n 
fa  a ra  a e 


k-  n n » k n " 
fe  a e re  e 


k n > k.  n n 
fa  a fe  a e 
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The  forward  reaction  rates  become  the  dominant  process  in  the  flow. 

For  cases  analyzed  in  Chapters  6 and  7 the  degree  of  ionization  for 
frozen  flow  is  larger  than  that  for  nonequilibrium  flow  at  fixed  x and  t|  values. 
As  x increases,  degree  of  ionization  increases  for  frozen  flow  and  decreases 
for  nonequilibrium  flow  at  fixed  ij  value. 

V ' 


APPENDIX  E 


EFFECTS  OF  BOUNDARY  LAYER  ON  SHOCK  WAVE  STRUCTURE 

The  analyses  in  Refs.  2k  and  25  on  shock- wave  structure  were  made 
with  the  assumption  that  flow  was  one-dimensional.  The  role  of  the  boundary- 
layer  growth  on  shock-vi:  ve  structure  was  not  considered.  Mirels'(Ref . 54)  has 
shown  that  the  flow  between  the  shock  and  contact  surface  in  an  actual  shock 
tube  is  nonuniform  due  to  the  side-wall  boundary  layer.  Recently,  Enomoto 
(Ref . 55)  has  studied  the  effects  of  the  boundary-layer  growth  on  the  ioniza- 
tion relaxation  in  argon,  based  on  MLrel's  boundary- layer  theory.  Enomoto 
(Ref.  55)  showed  that  the  temperature , density  and  pressure  increase  in  value 
with  distance  from  the  shock  front  due  to  boundary  layer  growth,  and  the 
ionization  relaxation  time  is  significantly  shortened  by  introducing  side-wall 
boundary- layer  effects. 

The  present  section  only  discusses  the  governing  equations  which  are 
applicable  to  quasi -one- dimensional  nonequilibrium  flow  behind  a shock  wave. 
Details  of  the  effects  of  the  side-wall  boundary  layer  on  the  shock-wave 

structure  is  under  study  and  will  be  presented  in  a forthcoming  UTIAS  Report. 

• 

The  specification  of  a problem  in  the  field  of  gasdynamics  requires 
the  flow  equations  (mass,  momentum  and  energy)  with  supplementary  information 
on  the  equation  of  state  of  the  gas.  The  governing  equations  for  flow  behind 
the  shock  wave  are  written  in  shock- fixed  coordinates  as 


I b <A  P u>  = 0 (=■) 

I |Z  <A  P u2>  = - i (E2) 

ih  p u H> = -s  (E3) 

i ^ <A  p u a>  = ma  ne  (E4) 

X £c  <A  P u he>  =<u  3T>+6el  + <E5> 


where  A is  the  effective  cross-sectional  tube  area,  and  <f>  = / f dA/A. 

By  using  the  equation  of  state,  p = pR(Ta  + QTe) , the  governing 
equations  dan  be  approximated  as  (Enomoto,  Ref.  55) 


E-l 


_ 2 %l  + Qinel 
<3x  ' 3 “ ne  u 


dl 

a 

dx 


u du 
C dx 


a 


dl 

e 

dx 


fe  du 

3u  dx 


Te  + 


5 AI  ) dx 


52  = 

dx 


du 

eud£ 


(e8) 

(E9) 

(E10) 


The  above  five  equations,  Eqs.  (E6  - E10) , are  required  for  the  five 
dependent  variables:  a,  u,  Ta,  Te  and  -i.  The  initial  conditions  at  x = 0, 
immediately  behind  the  frozen  shock,  are 


a « 0 


u 


u = 


T = 
a o 


M 


5 

TB 


M 


5M„ 


T„  cs  T or  T 
e a o 

P u 

p = -2-2  R T 
u a 


+ b] 


where  the  subscript  o denotes  quantities  evaluated  in  the  unionized  (a  = 0) 
upstream  gas.  The  length  required  to  reach  quasi-equilibrium  is  nearly 
identical  for  both  initial  values  of  Te. 

The  values  A and  dA/dx  can  be  obtained  from  the  solutions  of  side-wall 
boundary  layer  (Chapter  7) . 


It  has  been  shown  by  Glass  and  Liu  (Ref.  24)  that  the  radiation- energy 
loss  has  no  effect  on  the  relaxation  length.  The  condition  that  the  side-wall 
boundary  layer  effects  on  shock  structure  can  be  neglected  is  only  good  when 


dA  2p  R Tj.  A ^ 
dx  5p  dx 


(Ell) 


It  can  be  shown  that  the  side-wall  boundary-layer  effects  in  the 
UTIAS  10  cm  x 18  cm  Hypervelocity  Shock  Tube  at  shock  Mach  nuntoers  of  Mg  « 13  ~ l6  is 
not  small  since  Eq.  (Ell)  is  not  satisfied. 

A detailed  study  on  the  mutual  interactions  between  shock  structure  and 
shock-tube  sidewall  boundary  layer  flows  based  on  the  correct  set  of  effective 
quasi-one-dimensional  flow  equations  will  be  reported  by  K.  Takayama  and  W.  S.  Liu. 
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